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ABSTRACT 


A  new  generation  of  inexpensive  GPS  receivers  is  emerging, 
which  utilize  temperature  compensated  crystal  oscillators 
(TCXO)  as  clocks.  The  MX4200  has  been  chosen  as  a 
representative  sample  of  this  emerging  technology  receiver.  A 
detailed  discussion  on  the  receiver  architecture  and  its 
implications  in  the  measurement  process  is  given.  Experiments 
on  the  sensitivity  of  its  clock  to  environmental  effects,  as 
well  as  the  duration  of  the  transient  regime  are  described. 
One  limitation  to  the  use  of  these  receivers  in  geodetic 
applications  is  the  stability  of  its  clock.  The  noise  in  the 
receiver  clock  has  been  characterized  by  the  Allan  variance 
for  sampling  intervals  up  to  4000  seconds.  The  implication  of 
the  noise  levels  for  short  time  intervals  on  the  applicability 
of  the  receiver  to  high  precision  work  is  examined. 
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I.  IirrRODDCTIOM:  THE  SUBJECT  OF  THE  RESEARCH 

A.  BACK6R0U1ID 

Accurate  real  time  position  and  velocity  is  an  essential 
piece  of  knowledge  in  many  fields  of  science  and  engineering. 
Electronic  radionavigation  systems  (Decca,  Loran,  Omega,  etc.) 
have  been  successfully  used  in  the  last  decades.  Transit 
doppler  satellite  system  was  the  first  operational  spaceborne 
positioning  system  available  for  scientists  and  surveyors  in 
the  late  sixties. 

None  of  them  has  had  as  wide  an  application  as  the  Global 
Positioning  System  (GPS) .  First  conceived  as  a  military 
oriented  system,  its  impact  on  navigation  and  surveying  has 
been  remarkable.  Its  level  of  accuracy,  for  short  periods  of 
observation,  exceeds  any  other  available  system.  But  its 
success  relies  on  the  quality  of  its  main  component:  the 
clock. 

Different  kinds  of  clocks  have  been  used  for  this  purpose, 
mainly  quartz  oscillators,  rubidium  and  cesium  resonators. 
Among  other  things,  clocks  have  been  the  limiting  factor,  both 
in  price  and  size,  for  certain  applications.  Technology  in 
this  field,  as  well  as  in  integrated  circuitry,  has  improved 
dramatically  in  the  last  few  years.  Nevertheless,  crystal 
oscillators  are  very  sensitive  to  changes  in  temperature. 
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Usually,  manufacturers  place  the  oscillator  inside  an  oven  at 
a  fairly  high  temperature.  This  imposes  severe  restrictions 
regt  ding  power  supply  requirements. 

he  last  generation  of  surveying  receivers,  small  sized 
and  with  built-in  power  supply,  was  available  in  the  late 
eighties.  They  have  tracking  loops  that,  again,  require  a  high 
performance  quartz  oscillator.  They  are  expensive  devices, 
then,  although  capable  of  levels  of  accuracy  only  expected  for 
receivers  equipped  with  atomic  clocks  before. 

Great  improvements  have  occurred  in  the  design  of  a 
different  kind  of  receiver.  A  new  line  of  products  is  now 
emerging.  They  have  a  temperature-compensated  crystal 
oscillator,  not  enclosed  in  an  oven.  Both  price  and  power 
requirements  have  dropped  dramatically.  They  may  not  be  suited 
for  high  level  geodetic  tasks;  but  are  entirely  adequate  for 
certain  surveying  or  oceanographic  applications. 

Magnavox  released  its  last  GPS  receiver,  model  MX4200,  in 
January  1990.  It  weighs  less  than  a  kilogram,  costs  just  over 
$3,000  and  requires  less  than  3.5  watts,  D.C.  power. 
Experiments  on  its  performance  seemed  to  be  worthwile. 

B.  OBJECTIVE 

The  objective  of  this  research  is  the  characterization  of 
the  clock  performance  in  the  Magnavox  4200  GPS  receiver.  To 
this  end,  several  laboratory  tests  were  performed  to  obtain 
the  receiver  measurements  free  of  as  many  systematic  effects 
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as  possible.  An  analysis  of  the  stochastic  noise  was  then 
feasible. 

As  usual  with  brand  new  pieces  of  equipment,  a  great  deal 
of  research  had  to  be  devoted  to  try  to  understand  the 
receiver  idiosyncrasies,  which  turn  out  to  be  innovative  in 
many  respects.  As  emerging  technology,  it  was  subject  to 
software  and  firmware  modifications  to  upgrade  its 
performance.  All  the  efforts  have  been  done  to  update  the 
analysis  to  the  last  available  version. 

C.  RELATED  WORK 

Clynch  and  Coco  [Ref.  1]  have  investigated  the  short  term 
phase  noise  in  several  clocks,  from  GPS  measurements,  to 
determine  the  limits  of  reconstruction  of  the  number  of  whole 
cycles  lost  when  phase  lock  is  lost  and  regained.  Colquitt  and 
Anderle  [Ref.  2]  described  the  results  of  a  simulation  of  the 
effects  of  the  frequency  stability,  for  various  clocks,  on  the 
Doppler  Transit  measurements. 
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II.  AN  OVERVIEW  ON  THE  GLOBAL  POSITIONING  SYSTEM 


A  general  background  on  the  GPS  system  is  introduced  here 
to  give  the  reader  a  better  understanding  of  later 
explanations.  Most  of  it  should  be  credited  to  Wells  [Ref.  3]. 

A.  THE  GPS  CONCEPT 

The  GPS  is  based  on  two  elementary  principles: 

•  Given  the  position  of  a  satellite,  it  is  possible  to 
compute  the  satellite-receiver  distance  by  measuring  the 
time  it  takes  for  a  signal  sent  by  the  satellite  to  get  to 
the  receiver. 

•  Given  the  frequency  of  a  signal,  it  is  possible  to  compute 
the  radial  component  of  the  relative  velocity  by  measuring 
the  Doppler  shift  in  the  incoming  signal. 

This  will  provide,  then,  two  independent  measurements  per 

space  vehicle  (SV) .  If  three  SV's  are  being  observed,  it  is 

possible  to  solve  for  three  spatial  coordinates  and  three 

components  of  the  velocity  vector. 

Actually,  this  is  true  if  no  error  exists  in  the  receiver 

clock.  It  may  be  fairly  stable,  but  its  measurements  are 

scaled  by  the  speed  of  light,  which  is  a  huge  number.  That  is 

why  both  offset  and  drift  of  the  receiver  clock  are  computed 

as  a  part  of  the  navigation  solution.  This  is  achieved  by 

observing  either  four  SV's,  or  only  three,  with  an  added 

constraint:  the  antenna  elevation,  which  is  sometimes  fairly 
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well  known.  In  both  cases  a  deterministic  solution  can  be 
obtained: 

•  Three  SV^s,  range  and  Doppler  shift  each,  are  enough  to 
solve  for  x,  y,  x,  x,  y,  i. 

•  Four  SV's,  range  and  Doppler  shift  each,  provide  data  to 
solve  for  x,  y,  z,  x,  x,  y,  i,  f. 

If  more  than  four  satellites  are  available,  redundant 
observations  can  be  made.  It  leads  to  a  better  estimate,  via 
a  least  squares  solution. 

The  solution  is  usually  implemented  through  a  Kalman 
filter.  Current  measurements  update  a  solution  that  accounts 
for  its  recent  history.  This  allows  for  continuous  solution  in 
case  a  satellite  is  temporarily  lost  or  if  changes  in  the 
available  constellation  occur.  In  addition,  this  filters  out 
a  great  deal  of  the  noise  in  the  current  measurements. 

B.  THE  SIGNAL 

Three  independent  pieces  of  information  are  then  needed  to 
achieve  the  measurement: 

•  A  reference  signal  to  compare  with  the  incoming  Doppler 
shifted  signal. 

•  Something  to  recognize  when  the  signal  left  the  SV. 

•  The  position  and  dynamics  of  the  SV  at  that  precise  epoch. 
All  this  information  is  embedded  in  the  signal  the  SV 

sends:  it  consists  of  a  carrier  wave  of  a  known  frequency, 
which  is  modulated  by  'code'  and  message.  The  code  is  a 
repeating  trend  of  two  level  states  in  a  precise  pattern,  one 
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for  each  satellite.  It  begins  at  a  certain,  known  epoch  in  the 
SV  clock.  The  message  is  a  set  of  digitally  coded  pieces  of 
information  regarding  the  SV  position  and  status,  among 
others . 

Then  the  incoming  signal  can  be  expressed  as: 

S(t)=  C(t)  D(t)  exp(-i27r{f-fj)t)  (1) 

where  C(t)  represents  the  code,  D(t)  the  data  and  f-f^  is  the 
carrier  nominal  frequency,  shifted  by  a  Doppler  frequency  of 
fa- 

The  way  in  which  the  carrier  is  modulated  by  code  and 
message  is  not  pertinent  to  the  present  discussion.  But  all 
the  timing  processes,  as  well  as  the  signal  generation,  are 
driven  by  the  main  oscillator  in  the  SV.  Its  nominal 
frequency,  10.23  Mhz.  Different  scale  factors,  always 

integers,  provide  all  the  frequencies  and  time  tagging  used  in 
the  whole  process.  Thus,  an  error  detected  in  any  of  them 
leads  to  a  precise  knowledge  of  the  error  in  the  main 
oscillator,  and  in  turn,  in  any  of  the  frequency  and  time 
tagging  parameters. 

An  overview  of  all  three  follows: 

1.  Carrier 

Two  frequencies  are  transmitted  by  the  SV: 

•  LI  at  1575.42  Mhz  (154  f^) 

•  L2  at  1227.60  Mhz  (120  fg) 
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The  reasons  for  this  choice  are  discussed  in  Spilker 
[Ref.  4:p.  37].  Ionospheric  effects  are  minimum  over  this 
range  of  frequencies,  among  those  available  in  the 
radionavigation  bands.  At  the  same  time,  they  are  far  enough 
from  each  other  so  that  dual  frequency  measurements  can  be 
accurately  done,  as  discussed  later. 


Figure  l.-  Integrated  carrier  phase  generation 


The  Doppler  measurement  is  obtained  by  first  removing 
the  carrier  frequency  from  the  incoming  signal  (Figure  1) . 
Then  a  Doppler  tracking  loop  follows.  A  modeled  Doppler 
frequency  is  generated  by  the  tracking  loop  oscillator,  and 
subtracted  from  the  remaining  signal.  The  oscillator  frequency 
is  continuously  corrected  by  the  error  signal  after  the 
Doppler  frequency  removal.  This  modeled  frequency  is  then  the 
Doppler  frequency  measurement.  This  is  how  much  the  frequency 
of  the  incoming  signal  is  offset  from  the  nominal  carrier 
frequency,  due  to  the  dynamics  of  both  satellite  and  receiver. 
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If  Doppler  cycles  are  counted  since  tracking  was 
locked,  a  quantity  called  *  integrated  carrier  phase'  is 
obtained.  This  quantity,  times  the  carrier  wavelength,  is 
called  Phase  derived  range.  It  is  a  measurement  of  the  range 
to  the  satellite,  biased  by  the  range  at  the  time  when 
tracking  was  locked  on.  This  quantity  has  a  deep  impact  in 
ranging  for  precise  positioning.  Once  that  bias  can  be 
determined,  measurements  are  more  accurately  done  based  on  the 
phase.  A  one  per  cent  of  the  total  cycle  can  be  easily 
determined  on  a  phase  measurement.  The  wavelength  of  the 
carrier  is  about  0.19  meters;  as  a  consequence,  the  noise  in 
this  kind  of  measurement  is  in  the  millimeter  range. 

2 .  Coda 

This  refers  to  a  digitally  generated  pseudorandom 
noise  (PRN) .  This  is  a  long  sequence  of  transitions,  between 
two  states,  generated  at  a  certain  rate  and  repeating.  It  is 
known  as  a  Gold  code.  When  its  length  is  very  large,  it 
approximates  the  true  random  code  behavior.  This  is: 

•  It  has  an  energy  density  spectrum  comparable  to  a  random, 
white  noise.  That  is  where  its  name  comes  from. 

•  The  autocorrelation  function  shows  a  triangular  peak,  of 
unit  height  and  twice  the  chip  length  wide,  at  no  time 
delay.  It  is  1/N^  elsewhere,  N  being  the  niunber  of  chips 
integrated.  Here  chip  refers  to  the  interval  between  clock 
pulses. 

•  The  cross  correlation  between  different  codes  of  the  same 
nature  is  1/N^  for  any  time  delay.  They  behave  as  two 
uncorrelated  white  noises. 
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These  features  make  the  PRN  code  very  well  suited  for 
ranging  applications.  It  carries  only  two  pieces  of 
information:  SV  identification  and  transmit  time.  It  is 
therefore  possible  to  spread  the  signal  in  a  wide  bandwidth  at 
very  low  levels  of  power,  without  compromising  its 
reliability. 

Every  satellite  generates  two  different  codes: 

•  C/A  (clear/acquisition)  code,  1023  chips  long,  generated 
at  1.023  Mhz  and  self  repeating  every  millisecond.  The 
Standard  Positioning  Service  (SPS)  is  implemented  by  means 
of  this  code.  It  modulates  both  LI  and  L2  frequencies. 

•  P  (precision)  code,  267  days  long.  It  is  truncated  to  one 
week  length,  and  repeated  then.  It  is  generated  at  10.23 
Mhz,  so  its  chip  length  is  I/IO  that  of  the  C/A  code.  This 
allows  levels  of  accuracy  an  order  of  magnitude  better 
than  using  C/A  code.  The  Precise  Positioning  Service  (PPS) 
is  implemented  through  this  code.  It  modulates  only  the  LI 
frequency. 

Implicit  in  the  code  is  the  transmit  time.  The  code 
generation  is  driven  by  the  main  oscillator  in  the  satellite. 
The  oscillator  cycles  count  is  in  fact  the  realization  of  the 
clock.  Thus,  each  transition  corresponds  to  a  precise  time  in 
the  satellite  clock. 

The  measurement  is  done  by  means  of  a  code  tracking 
loop.  A  replica  of  the  corresponding  code  is  generated  in  the 
receiver.  It  is  been  delayed  until  maximum  output  from  a 
correlator  is  obtained.  Once  this  process  is  locked  on,  the 
code  stays  tracked  by  feeding  back  the  corrections,  due  to 
misalignment  between  code  and  replica,  to  the  existing  delay. 
This  delay  is  then  the  essence  of  the  measurement,  which  is 
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requested  at  a  given  epoch  time  of  measurement.  This  delay, 
times  the  speed  of  light,  is  called  pseudoranoe . 

3 .  Data 

Generated  at  50  bps  chip  rate,  there  is  a  data  signal 
that  takes  12.5  minutes  to  be  completed.  Within  this  period, 
information  regarding  the  tracked  satellite,  termed  ephemeris . 
is  repeated  every  30  seconds.  The  ephemeris  data  contains  the 
orbital  parameters  of  the  SV.  It  allows  the  user  to  compute 
the  SV  position  and  velocity.  Corrections  to  the  received 
signal,  due  to  misalignment  between  the  SV  clock  and  the 
absolute  time,  are  contained  in  the  ephemeris,  too.  Thus,  data 
in  the  message  allows  the  user  to  know,  at  the  best  real  time 
accuracy  level,  the  position  of  the  satellite  and  the  error  in 
its  clock.  The  high  rate  at  which  these  data  are  transmitted 
allows  the  reciver  the  acquisition  of  a  satellite  in  a  fairly 
short  time. 

The  rest  of  the  information,  termed  almanac,  is 
related  to  the  whole  constellation.  It  facilitates  the 
searching  of  any  other  satellite.  Status  of  the  whole  system 
is  broadcast  here,  too,  as  well  as  many  other  special 
messages. 

C.  SYSTEM  SEOMEMTS 

It  is  customary  to  divide  the  system  in  three  main 
components:  control,  space  and  user  segments. 
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1.  Control  sognent 

This  consists  of  five  monitoring  stations  around  the 
globe,  which  continuously  track  all  the  satellites  in  view.  By 
fitting  the  recently  measured  ranges  to  an  orbital  model,  the 
Operation  Control  Center  (OCS)  adjusts  the  parameters  of  that 
orbit.  That  data  is  injected  to  the  SV's,  as  a  new  ephemeris 
set,  daily.  These  ephemeris  sets  are  considered  valid  for  two 
hours  before  and  after  their  time  of  epoch.  They  give  the 
position  of  the  SV  within  1  part  per  million.  Almanacs  are 
recomputed  and  injected  to  the  whole  constellation  on  a  weekly 
basis. 

Time  keeping  is  the  second  major  issue  in  the  control 
segment.  Every  control  station  is  equipped  with  cesium  atomic 
clocks;  even  so,  they  differ  relative  to  each  other. 
Currently,  five  cesium  clocks  at  the  OCS  serve  as  the  master 
clock.  They  are  adjusted  every  two  days  by  comparison  with  the 
standard  at  the  US  Naval  Observatory.  Once  the  ground  stations 
have  an  absolute  time,  it  is  possible  to  solve  for  the  SV 
clock  error.  A  second  order  polynomial  model  of  its  clock  is 
included  as  a  part  in  the  ephemeris  set. 

It  is  worthwhile  to  say  something  about  the  GPS  time 
system.  It  was  set  to  zero  at  OOh  00m  00s  UTC,  January  6th 
1980,  and  is  based  on  weeks  of  604800  seconds  of  UTC  time.  The 
time  is  then  expressed  in  weeks  (WN)  and  seconds  within  the 
week  (SOW) .  Because  of  the  irregular  rate  of  rotation  of  the 
Earth,  the  time  community  corrects  the  deviation  between  UTC 
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and  true  rotation  by  adding/subtracting  a  leap  second  once  in 
a  while.  This  is  not  true  for  GPS  time.  That  is  why  there  is 
a  7  second  difference  between  GPS  and  UTC  times  as  of 
September  1991. 

2 .  Space  segment 

This  consists  of  the  satellites  themselves.  As  stated, 
the  SV's  broadcast  the  almanac  and  ephemeris  data,  in  addition 
to  generating  C/A  and  P  codes.  All  electronics  are  driven  by 
a  main  atomic  oscillator  (they  have  spares,  too)  whose  error 
is  reported  in  the  ephemeris. 

The  orbital  geometry  has  been  designed  to  provide  3-D 
continuous  coverage  worldwide:  six  orbital  planes,  at  55° 
inclination,  with  four  SV's  each.  The  orbit  is  almost  circular 
and  the  height  above  the  Earth's  surface  is  about  20,000  Km. 
This  makes  the  orbital  period  12  hours.  At  such  a  height, 
about  5  hours  out  of  12,  a  satellite  is  above  the  observer's 
horizon. 

As  of  September  1991,  only  17  SV's  are  available. 
This,  indeed,  provides  2-D  continuous  coverage  worldwide.  A  3- 
D  coverage  is  expected  by  the  end  of  1991. 

3 .  User  segment 

This  represents  the  user's  end  of  the  system:  the 
receiver.  A  great  variety  of  receivers  have  been  designed,  and 
they  differ  in  many  respects: 

•  Single  (LI)  versus  dual  frequency  (LI  and  L2) . 


12 


•  C/A  versus  C/A+P  codes  tracking. 

•  Main  oscillator:  quartz,  rubidium,  cesium  and 

combinations. 

•  Single  tracking  channel  versus  multiple  channels. 

•  Continuous  tracking  versus  sequential  and  multiplexing. 

All  receivers  have  in  common  the  existence  of  an  error 
in  their  clocks,  which  has  to  be  computed  as  a  part  in  the 
navigation  solution,  regardless  how  small  it  is.  The  final 
application  usually  determines  its  performance:  it  is  related 
to  the  desired  level  of  accuracy,  as  well  as  other  constraints 
in  size,  power  requirements  and  price. 

D.  BIASES  AMD  ERRORS  IM  THE  KEASUREMEMT 

As  usual  in  the  scientific  literature,  a  bias  is  an  error 
in  a  measurement  that  can  be  identifiable  and  somehow  modeled 
or  computed.  It  may  or  may  not  be  removed  from  the 
measurement.  A  true  error,  on  the  other  hand,  can  have  an 
identifiable  source,  but  cannot  be  modeled.  It  properties  can 
be  estimated.  Statistical  variables  are  used  to  describe  its 
effect  on  the  measurement  accuracy. 

Biases  and  errors  must  be  properly  treated  in  making  a  GPS 
measurement.  They  can  be  classified  into  four  main  groups, 
depending  on  their  source:  due  to  the  SV,  to  the  propagation 
medium,  to  the  receiver  and  observation  errors. 
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1.  Biases  due  to  the  satellite 

There  are  two  sources  of  error  here: 

a.  Bias  in  the  SV  clock 

A  model  of  it  is  contained  in  the  ephemeris 
message,  so  a  great  deal  of  it  can  be  removed  as  a  bias.  This 
is  effectively  done,  in  software,  by  all  the  receivers.  Some 
small  errors  (up  to  10  m)  remain. 

Jb.  Bias  in  the  SV  position 

The  real  time  solution  is  produced  by  using  the 
orbital  parameters  in  the  broadcast  ephemeris.  Those 
parameters  correspond  to  a  projected  orbit,  as  computed  from 
previous  range  data  to  the  OCS  stations.  Thus,  it  is  based  on 
an  extrapolation  technigue.  How  well  those  parameters  will 
describe  the  true  orbit  determines  the  size  of  this  bias.  It 
can  be  estimated  about  5m,  sometimes  up  to  10. 

A  smaller  bias  can  be  obtained  using  what  is  called 
'precise  ephemeris'.  This  set  of  orbital  parameters  is 
computed  after  tracking  a  SV  for  eight  days.  It  is  available 
for  a  postprocessing,  and  is  based  on  an  interpolation 
technique.  Precise  ephemeris  provides  a  more  accurate  SV  clock 
model,  too.  The  Defense  Mapping  Agency  (DMA)  and  the  National 
Geodetic  Survey  (NGS)  have  their  own  tracking  station 
networks.  Both  agencies  produce  precise  ephemeris,  but  its 
distribution  is  very  limited. 
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2.  Biases  due  to  the  propagation  medium 

There  are  two  sources  of  error  here: 


a .  Ionospheric  delays 

The  ionosphere  is  a  dispersive  medium  for  an 
electromagnetic  wave.  The  ultraviolet  light  ionizes  the 
molecules  in  the  ionosphere,  releasing  free  electrons.  These 
electrons  have  the  effect  of  changing  both  the  phase  and  group 
velocity  of  the  wave.  For  a  differential  element  of  time,  dt, 
the  local  element  of  measured  range  will  be 

dfi=c  dt  ( 2 ) 

because  a  velocity  of  c  is  assumed.  For  the  same  element  of 
time,  the  true  range  is 

dr=v  dt  (3) 


because  v  is  the  actual  velocity  of  propagation. 

The  ratio  between  those  elements  is  defined  as  the 
index  of  refraction.  It  is,  in  fact,  the  ratio  between  assumed 
and  true  velocities: 


do  c 
n  =  ^®_ 
os  V 


(4) 


The  measured  range,  then,  will  be 
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P  S 


p  = 


ds 


(5) 


Pseudorange  is  a  measurement  of  a  time  delay.  Thus, 
it  is  based  on  the  propagation  of  the  wave  energy.  The  energy 
propagates  at  the  group  velocity.  On  the  other  hand,  the  phase 
of  the  signal  travels  at  the  phase  velocity.  Local  group  and 
phase  indices  of  refraction  can  be  expressed,  as  a  high 
frequency  approximation,  as  [Ref.  3:p.  9-2]; 


N 


(6) 


n  =  — =  (l-aiL) 


(7) 


where  v^,  are  the  group  and  phase  propagation  velocity  and 
n^,  rip  group  and  phase  indices  of  refraction,  a  is  a  certain 
constant,  f  the  frequency  and  N  is  the  ionospheric  electron 
density  (electrons/m’)  . 

The  measured  pseudorange,  thus,  will  be 

p=r+-f-fN  ds=r*a^  (®) 


whereas  the  measured  phase  range  is 
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(9) 


N,  refers  to  integrated  electron  count  along  the  path. 

Estimated  quantities  for  these  variables  are: 

•  a  about  40  m^  s'^. 

•  N  about  lO'^  m'^  (electrons  per  cubic  meter) . 

•  f  about  1.5  X  10*  Hz  for  the  LI  carrier. 

N,  has  a  maximum  in  daytime,  and  a  minimum  during 
the  nighttime.  It  depends  also  on  the  elevation  angle:  the 
lower  the  SV,  the  larger  the  integration  path  through  the 
ionosphere . 

Calling  d,  =  (a  the  ionospheric  delay,  the 

measured  pseudorange  is  long  by  d,,  and  the  measured  phase 
derived  range  is  short  by  the  same  amount.  Spilker  [Ref.  4:p. 
35]  estimates  that  this  delay,  for  vertical  incidence,  is  10 
ns  (3m)  during  nighttime  and  up  to  50  ns  (15m)  in  daytime. 
About  three  times  these  quantities  occur  for  elevation  angles 
of  10  degrees. 

This  quantity  is  hard  to  model.  But  it  can  be 
easily  estimated  from  dual  frequency  observations  by  taking 
the  difference  between  phase  range  measurements: 
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then 


^L2 


(11) 


b .  Tropospheric  delays : 

The  troposphere,  on  the  other  hand,  is  not  a 
dispersive  mediiun  for  the  GPS  range  of  frequencies.  It  affects 
both  phase  and  range  measurements  in  the  same  way.  It  is 
caused  by  water  vapor  molecules,  as  well  as  other 
constituents.  Again,  its  effect  is  minimum  at  highest 
elevation.  The  dry  component  is  about  90%  of  the  total,  and 
can  be  accurately  estimated  from  surface  atmospheric  pressure 
measurements.  The  wet  component  is  not  so  accurately  known, 
because  water  vapor  content  is  extremely  variable  with  height. 
Its  effect  is  much  smaller,  though.  There  are  several  models 
currently  available  to  correct  for  tropospheric  delays.  Most 
of  the  geodetic  receivers  implement  an  algorithm  to  correct 
it,  based  on  meteorological  data  input. 

At  sea  level  and  vertical  incidence,  a  rough 
approximation  of  the  total  tropospheric  range  error  is  2.3 
meters  [Ref.  4:p.  36].  It  can  reach  up  to  20  meters  at  a  10 
degree  elevation  angle. 

3.  Biases  and  errors  due  to  the  receiver 
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a.  Receiver  clock 


This  is  the  main  source  of  errors  in  the 
measurement.  The  time  difference  between  the  receiver  clock 
and  GPS  time  produces  a  range  error  which  is  that  difference 
multiplied  by  the  speed  of  light.  The  frequency  difference 
between  the  receiver  clock  (oscillator)  and  the  nominal 
carrier  frequency  produces  a  velocity  error  of  that  difference 
multiplied  by  the  carrier  wavelength. 

Most  receivers,  as  stated,  solve  for  time  and 
frequency  biases  as  a  part  of  the  navigation  solution.  These 
errors  can  be  used  then  to  correct  the  actual  measurements  by 
their  range  and  velocity  equivalents.  How  well  this  adjustment 
is  achieved  depends  strongly  on  the  stability  of  the  clock,  as 
well  as  on  the  quality  of  the  navigation  filter. 
b.  Carrier  phase  ambiguity 

Most  of  the  high  performances  receivers  use  the 
integrated  carrier  phase  to  improve  the  range  measurement.  The 
phase  measurements  are  more  accurate  than  the  delay 
measurements . 

Wells  [Ref.  3:p.  9.14]  estimates  the  resolution  of 
the  measurement  to  be  about  one  percent  of  the  signal 
wavelength.  The  C/A  code  chip  length  is  about  300  meters,  and 
the  P  code  30  meters.  This  leads  to  residual  errors  on  the 
order  of  3  and  0.3  meters  respectively.  LI  frequency,  on  the 
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other  hand,  has  a  wavelength  about  20  cm,  so  residual  errors 
could  be  about  2  mm,  within  a  given  cycle. 


Integrated  carrier  phase  is  a  biased  measurement, 
though.  The  true  range  to  the  SV  at  time  of  lock-on  is 
unknown:  there  is  a  certain  integer  number  of  cycles  that  has 
to  be  resolved.  There  are  many  different  approaches  to  this 
problem.  If  the  estimate  of  that  bias  is  off  by  N  cycles,  this 
goes  as  a  whole  to  the  measurement,  as  an  uncorrected  bias. 
This  bias  amounts  N  times  the  carrier  wavelength. 

The  same  problem  is  found  in  a  situation  when 
tracking  of  a  given  satellite  is  lost.  When  tracking  is 
resumed,  an  integer  number  of  wavelengths  have  been  lost  from 
their  count.  This  is  known  as  'cycle  slips'.  Different 
techniques  have  been  developed  to  deal  with  this  problem.  If 
that  count  is  not  properly  retrieved,  the  solution  remains 
biased. 

c.  Electronics,  processing  and  computational  noise 
There  are  a  number  of  noise  sources  along  the 
signal  path  in  the  receiver:  antenna,  pre-amp,  demodulator..., 
all  of  which  generate  thermal  noise.  The  resolution  of  the 
measurement  is  subject  to  hardware  limitations;  these  lead  to 
residual  processing  noise,  too.  Finally,  there  is  noise  in  the 
computation  process  in  the  navigation  filter:  round-off 
residuals,  etc. 
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Electronics  and  computational  noise  are  considered 
to  be  at  a  much  lower  level  than  any  other  error  here 
mentioned.  Processing  noise  in  the  pseudorange  measurement  can 
be  fairly  high,  though.  Because  of  that,  precise  positioning 
is  done  based  on  phase  range  measurements. 

4 .  Observation  errors 

Here  is  a  group  of  true  errors  which  cannot  be 
removed,  although  some  of  them  can  be  easily  detected. 

a.  Multipath 

Metallic  surfaces  act  as  a  mirror  for  the  incoming 
RF  signal.  What  the  receiver  gets  is  a  summation  of  direct  and 
reflected  signal.  This  has  the  effect  of  shifting  the  phase  of 
the  signal  by  a  certain  amount,  depending  on  the  relative 
geometry  (SV-ref lector-antenna) .  So,  it  shows  an  interference 
pattern  as  the  geometry  changes  and  is  easy  to  recognize. 

b.  Antenna  phase  center 

The  radiation  pattern  of  the  antenna  is  not 
radially  symmetric,  as  it  should  be.  In  fact,  each  antenna 
shows  a  different  radiation  pattern.  So,  the  phase  center  is 
dependent  upon  the  azimuth  and  elevation  of  the  tracked  SV. 

It  is  very  hard  to  know  the  true  radiation  pattern. 
How  well  it  is  known  determines  the  size  of  the  residual 
error.  It  is  usually  at  the  millimeter  level. 
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E.  SELECTIVE  AVAILABILITY  AND  LEVELS  OF  ACCURACY 

As  mentioned  above,  two  different  services  are  offered  by 
the  system:  Precise  Positioning  Service  (PPS)  and  Standard 
Positioning  Service  (SPS) .  Both  of  them  are  available  for 
military  and  certain  civilian  users.  Most  civilian  users, 
however,  must  use  the  SPS.  They  differ  in  the  levels  of 
accuracy  because  of  the  size  of  their  code  chips:  300  meters 
for  C/A  (SPS) ,  30  meters  for  P  code  (PPS) . 

The  system  was  first  designed  for  military  applications. 
US  DoD  reserves  the  ability  to  degrade  both  signal  and 
message,  so  the  quality  of  the  solution  gets  worse.  This  is 
known  as  the  'Selective  Availability'  (SA)  option.  This  might 
be  done  by  broadcasting  inaccurate  ephemeris  and  by 
intentionally  adding  noise  to  the  codes.  In  any  event,  the 
quality  of  the  signals  guarantees  the  following  levels  of 
accuracy : 

•  100  meters,  2-D  rms  for  SPS. 

•  16  meters,  3-D,  1  a,  for  PPS. 

In  an  SA  situation,  the  above  accuracy  of  PPS  requires 
special  devices  to  de-crypt  the  degraded  signals.  Thus,  it  is 
useless  for  civilian  users. 

Related  to  the  levels  of  accuracy  is  the  concept  of 
Geometrical  Dilution  of  Precision  (GDOP) .  It  represents  how 
much  the  precision  in  the  solution  is  degraded  due  to 
geometrical  considerations.  In  other  words,  how  degraded  the 
actual  solution  can  be,  as  compared  to  an  idealistic  situation 
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where  all  tracked  satellites  are  in  an  optimal  spatial 
distribution. 

The  variation  in  the  measured  range  to  the  j-th  SV,  due  to 
changes  in  the  receiver  position  and  clock  offset  can  be 
expressed  as: 


(12) 


where  the  partial  derivatives  are  only  determined  by 
geometrical  considerations.  For  measurements  to  n  SVs,  a 
system  of  equations  is  foznned: 

dp=/f  dx  (13) 

where  each  element  of  the  design  matrix  H  is 

(14) 


In  solving  the  system,  the  solution  is  found  as: 

dx  =  dp  (15) 

is  called  the  covariance  matrix.  Its  diagonal  elements 
represent  the  variance  in  the  corresponding  component  of  the 
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solution.  Again,  only  geometrical  considerations  affect  their 
size. 

GDOP  is  defined  as  the  square  root  of  the  trace  of  the 
covariance  matrix.  For  practical  considerations,  GDOP  is 
usually  'projected'  into  3  spatial  coordinates  and  time 
components.  Those  components  are  termed  'Horizontal  Dilution 
of  Precision'  (HDOP) ,  'Vertical  Dilution  of  Precision'  (VDOP) , 
and  'Time  Dilution  of  Precision'  (TDOP) .  This  is,  taking  only 
the  first  two,  the  third  and  the  fourth  terms  of  the  diagonal 
elements  of  respectively. 

GDOP  can  get  a  minimum  value  of  one.  The  larger  GDOP  is, 
the  larger  the  variance  in  the  solution  components.  This  is, 
the  accuracy  in  the  solution  gets  worse,  just  due  to  the 
relative  geometry  SV's-receiver. 

F.  DIFFEREMTIAL  POSITIONING 

When  the  GPS  system  was  first  designed,  it  was  not 
envisioned  that  there  would  be,  in  fact,  several  ways  to 
overcome  the  limitation  in  the  accuracy  levels.  It  was 
conceived  that  a  user  would  compute  a  real  time  solution, 
based  on  pseudorange  measurements  to  determine  the  position. 

There  are  many  receivers  that  include  phase  measurements 
in  the  computation  of  the  position  nowadays.  In  addition  to 
that,  several  techniques  have  been  found  to  improve  the 
solution  in  a  post  processing.  Still,  there  is  a  technic[ue 
that  allows  the  common  user  to  improve  a  real  time  solution: 
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differential  positioning.  It  is  based  upon  the  assumption  that 
most  of  the  errors  in  the  measurement,  external  to  the 
receiver  itself,  are  common  to  two  receivers  close  enough  to 
one  another.  A  reference  receiver  can  be  placed  on  a  well 
surveyed  point.  By  comparing  its  solution  with  the  true 
position,  corrections  either  to  measurements,  or  to 
coordinates,  can  be  computed.  They  are  then  sent,  real  time, 
to  the  second  receiver,  which  automatically  corrects  its  own 
solution.  This  is  a  very  powerful  technique,  of  common  use  in 
the  surveyor's  community.  It  is  being  used  for  navigation 
purposes ,  too . 

Differential  techniques  are  used  too  in  testing  receivers. 
If  the  position  of  two  receivers  is  known,  the  difference 
between  simultaneous  measurements  eliminates  most  of  the 
biases  and  errors  that  are  external  to  the  receivers.  The 
distance  between  receivers  has  to  be  small  enough  to  assume 
that  ionospheric  and  tropospheric  delays  are  the  same  for  both 
receivers.  Those  differences  provide  an  estimate  of  the 
receiver  performance. 
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III.  THE  MXGNAVOZ  4200  RECEIVER 


A.  GENERAL  CHARACTERISTICS 

The  MX  4200  has  been  designed  for  navigation  purposes, 
although  its  use  for  precise  positioning  is  the  object  of 
intense  research  [Ref.  5].  It  is  a  single  frequency  receiver, 
restricted  to  SPS  users.  It  can  continuously  track  the  C/A 
code  and  LI  frequency  of  up  to  six  satellites.  The  choice  of 
satellites  to  be  tracked  is  automatically  done  by  analyzing 
the  GDOP  each  possible  combination  produces. 

This  receiver  contains  software  for  real  time  navigation. 
The  navigator  consists  on  an  8  state  modified  Kalman  filter. 
It  provides  a  solution  in  three  spatial  coordinates,  three 
components  of  the  velocity  vector,  as  well  as  clock  offset  and 
drift  (frequency  offset) . 

The  receiver  clock  is  driven  by  a  temperature  compensated 
crystal  oscillator  (TCXO) .  This  kind  of  oscillator  is  very 
inexpensive  and  low  powered,  compared  to  oven-enclosed  crystal 
oscillators  and  more  sophisticated  clocks.  This  does  not 
degrade  the  real  time  solution.  The  navigation  filter  is  well 
designed,  and  clock  offset  and  drift  are  accurately  estimated, 
and  position  and  velocity  measurements  corrected  for  them.  But 
the  clock  behavior  can  be  a  major  issue  for  high  quality 
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measurements.  In  this  situation,  phase  derived  range,  rather 
than  pseudorange,  is  used  to  compute  the  solution  in  post¬ 
processing.  This  is  precisely  the  area  of  the  research  that 
has  been  done,  and  reported  here. 

The  receiver  also  can  work  in  a  differential  mode. 
Corrections  to  measurements,  rather  than  to  solutions,  can  be 
input  to  the  receiver,  following  the  RTCM-104  standard  format 
[Ref.  6]. 

This  receiver  was  first  conceived  as  a  GPS  upgrade  to  the 
Magnavox  4102  Transit  Navigator,  i.e.,  as  a  GPS  sensor  for  an 
existing  piece  of  hardware  [Ref.  7:p.  1].  Neither  control  nor 
data  logging  can  be  achieved  without  an  external  device.  But 
the  receiver  works  on  its  own,  automatically  providing 
solution,  as  well  as  raw  measurement  output.  The  next  section 
will  be  devoted  to  its  interfaces. 

The  above  mentioned  characteristics,  as  well  as  a  very 
high  level  of  integrated  circuitry,  allowed  the  manufacturer 
to  build  the  MX4200  at  a  very  low  cost  (just  over  $3000)  , 
small  sized  (18x14x4.5  cm,  0.9  Kg)  and  with  very  low  power 
requirements  (10  to  32  Vdc,  5  watts  max.)  [Ref.  7,p.  1]. 

The  antenna  is  a  quadrafilar  helix  in  an  all  weather  case 
(19  cm  high,  8.3  cm  diameter)  with  no  ground  plate.  For 
navigation  applications,  this  is  a  desirable  characteristic, 
to  avoid  loss  of  satellites  when  roll/pitch  occurs.  But  it  has 
its  trade-off  for  ground  applications,  because  it  is  more 
vulnerable  to  multipath  effects. 
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B.  IIITERFXCB8 


Four  different  interfaces  have  been  established:  Control 
port,  Raw  Data  port.  Equipment  port  and  RTCM  port.  In 
addition,  a  timing  pulse  can  be  output  once  a  second.  It 
allows  external  devices  to  reference  their  measurements  to  GPS 
or  UTC  time. 

1 .  Control  port 

This  is  an  I/O  interface  to  allow  the  control  of  the 
receiver  functions,  from  an  external  device.  The  navigation 
solution  is  output  through  this  port  as  well.  A  piece  of 
software  is  provided  by  the  manufacturer  to  perform  this 
operation  from  an  external  IBM-PC  computer,  connected  to  the 
receiver  through  the  control  port.  The  external  computer  then 
becomes  the  Control  and  Display  Unit  (CDU) .  Under  this 
configuration,  several  options  in  the  receiver  are  accessible 
to  the  user: 

•  Initialization:  input  of  approximate  coordinates  to 
facilitate  the  lock-on  process. 

•  Rejecting  undesired  satellites. 

•  Time  frame  selection  (GPS,  UTC) 

•  Navigation  mode:  2-D,  3-D,  altitude  aiding. 

•  Elevation  mask  for  tracked  satellites. 

•  Time  recovery  option. 

•  Requesting  ephemeris  and  almanac  data  in  the  GPS  message. 

•  Selection  of  desired  records  in  the  raw  data  stream,  to  be 
sent  to  the  raw  data  port. 
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•  Enabling/disabling  differential  corrections. 

•  Changing  the  dynamic  response  of  the  navigation  filter. 

Also,  a  continuous  display  of  the  status  of  the 
receiver,  tracked  satellites  and  navigation  solution  is 
available.  Current  navigation  solution,  differential 
corrections,  almanac  and  ephemeris  data,  can  be  also  recorded 
in  the  CDU. 

2.  Raw  data  port 

This  interface  allows  the  user  to  access  the  raw  data 
measurement,  as  well  as  navigation  solution  and  time  recovery 
data.  It  can  be  connected  to  an  external  IBM-PC  computer  to 
log  those  data.  Almanac  and  Ephemeris,  as  well  as  several 
status  and  debugging  messages,  are  output  through  this  port, 
too.  A  description  of  their  contents  and  meaning  is  given  in 
the  next  chapter. 

3 .  Equipment  port 

Allows  the  receiver  to  communicate  with  external 
navigation  devices,  following  the  NHEA-0813  Standard  for 
interfacing  marine  electronics.  Position,  velocity  and  time 
can  be  sent  to  those  devices  through  this  port. 

4 .  RTCM  port 

Inputs  the  corrections  to  the  measured  pseudoranges  as 
received  through  a  data  link,  following  the  standard  RTCM 
format.  Corrections,  then,  can  be  applied,  real  time,  to  the 
measured  pseudorange,  in  a  differential  navigation  mode. 
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C.  NAVIGATION  MOOES 


The  navigation  filter  can  compute  the  real  time  solution 
in  different  modes,  as  a  function  of  the  number  of  available 
SV's.  The  computation  of  the  components  of  the  state  vector  in 
the  Kalman  filter  is  prioritized  as  follows: 

If  4  SV's  are  being  tracked,  measurements  on  their  codes 
and  carrier  phase  provide  information  to  get  a  deterministic 
solution  for  the  entire  state  vector.  In  fact,  a  deterministic 
solution  is  only  obtained  at  initialization,  if  four  SV's  are 
available.  The  Kalman  filter  is  a  Least  Squares  Adjustment 
algorithm  that  uses  the  recent  history  of  the  solution  as 
extra  information  to  get  the  current  solution.  So,  even  if 
only  four  SV's  are  available,  the  solution  is  not  truly 
deterministic.  More  than  4  SV's  add  redundant  observations,  so 
the  reliability  of  the  solution  increases. 

The  user  is  allowed  to  fix  the  height  of  the  antenna, 
which  is  fairly  well  known  in  a  sea  environment.  Then,  three 
SV's  is  enough  to  get  a  deterministic  solution.  Height  and 
upwards  velocity  are  not  solved  for,  i.e.,  only  a  2-D  position 
and  velocity  solution  is  obtainable.  Transition  from  3-D  to  2- 
D  solutions  is  automatically  done  by  default,  if  the  receiver 
was  set  to  'auto'  altitude  aiding  mode. 

When  the  measurements  are  so  few  that  the  system  is  under 
determined,  the  Kalman  filter  keeps  a  'dead  reckoning'  (DR) 
navigation  mode. 
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The  selection  of  the  satellites  to  be  used  in  the  solution 
is  based  upon  limiting  values  of  elevation  mask,  VDOP  and  HDOP 
the  user  enters  as  a  part  xn  the  initialization  process.  In 
addition,  the  user  can  modify  the  dynamic  response  of  the 
navigator.  The  navigator  then  filters  the  measurements  as  a 
function  of  limiting  values  in  horizontal  and  vertical 
accelerations . 

D.  TIME  RECOVERY 

This  option  directs  the  receiver  to  take  the  measurements 
at  nearly  even  seconds  of  GPS  time.  This  is  done  by  adjusting 
the  local  physical  time  frame  by  the  results  of  the  Kalman 
solution.  In  addition,  a  timing  pulse  is  output  once  a  second. 

Measurements  are  taken  at  certain  epochs  in  the  receiver 
time  frame,  i.e.,  in  the  time  as  kept  by  the  receiver  clock, 
instead  of  the  GPS  time  frame.  This  is  known  as  a  'user  time 
receiver'  [Ref.  7:p.  3].  When  time  recovery  (TR)  is  disabled, 
measurements  are  taken  at  even  seconds  in  the  receiver  time 
frame.  When  enabled,  they  are  taken  as  close  to  an  even  GPS 
second,  as  the  receiver  is  able  to  determine,  based  in  the 
previous  solution.  This  makes  it  possible,  not  only  to  have 
the  timing  pulse  referenced  to  GPS  time,  but  also  to  know  the 
instantaneous  clock  offset  and  drift  at  that  precise  time.  The 
navigation  solution  corresponds  to  the  time  when  the  pulse  was 
output,  so  the  receiver  reports  what  was  the  status  of  its 
clock  at  the  very  instant  of  the  pulse  just  sent.  At  the  same 
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tine,  the  last  solution  is  used  to  compute  when,  in  the 
receiver  time  frame,  the  next  measurements  must  be  taken. 

B.  THE  MEASDREMEMT  PROCESS;  TRACKING  LOOPS 

The  receiver  operation  is  based  on  two  major  components: 
tracking  loops  and  the  navigation  filter.  Each  one  runs  on  its 
own,  with  an  interface  in  between  which  sends  the  measurements 
from  the  tracking  loop  to  the  navigator,  whenever  the  latter 
requests  them. 

1 .  Tracking  loops 

Most  of  the  information  here  contained  should  be 
credited  to  Keegan  [Ref.  7]  and  has  been  experimentally 
verified,  as  shown  in  next  chapter.  A  simplified  discussion 
will  be  presented  here. 

The  uppermost  part  of  Figure  2  shows  a  schematic  view 
of  the  tracking  loops.  At  a  first  stage,  carrier  frequency  is 
removed  from  the  incoming  signal.  It  is  generated  by  the  main 
oscillator,  at  the  nominal  LI  frequency.  The  frequency  error 
in  the  receiver  clock  is  then  transferred  to  the  Doppler 
frequency  that  is  left. 

A  replica  of  the  tracked  code  is  generated.  The  code 
is  then  removed  by  multiplying  that  replica,  properly  delayed, 
by  the  signal.  If  that  delay  was  well  computed,  the  cross¬ 
correlation  between  both  signals  is  one,  so  only  a  continuous 
wave  of  frequency,  modulated  by  the  data,  is  left. 
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TRACKING  LOCP 


OUTPUT 


Doppler  frequency  is  then  removed  by  multiplying  the 
signal  by  a  continuous  wave  of  fj  frequency.  Only  the  data  is 
left,  plus  an  additional  error  signal  due  to  the  mismodelling 
of  time  delay  and  Doppler  frequency.  By  integrating  this  error 
signal  in  the  correlator,  corrections  to  time  delay  and 
Doppler  frequency  are  computed.  Then,  they  are  fed  back  to  the 
code  generator  and  Doppler  frequency  oscillator.  This 
oscillator  also  drives  the  code  generator,  so  the  signal  is 
reconstructed  at  the  proper  chip  rate,  according  to  the 
Doppler  shift  in  the  incoming  signal. 

In  addition,  Doppler  frequency  is  integrated  once 
tracking  was  locked.  Scaled  by  the  carrier  wavelength,  it 
forms  the  current  'phase  range'  measurement.  This  integrated 
phase  is  then  applied  to  the  time  delay.  So,  the  code  replica 
is  generated  with  a  time  delay,  relative  to 
transmit  time  t^,  rather  than  t-t,,  as  usual  in  other 
receivers. 

In  this  way,  all  the  relative  dynamics  SV-receiver  has 
been  removed  from  the  code.  Tracking,  thus,  can  be  done  in  a 
much  narrower  bandwidth,  increasing  the  reliability  of  the 
measurement.  It  is  the  difference  between  pseudorange  and 
phase  range  that  the  code  loop  is  tracking.  Its  small  dynamics 
allows  for  filtering  the  code  signal  (which  is  much  noisier 
than  Doppler)  with  a  time  constant  quite  large.  This  kind  of 
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cycles.  The  next  measurement  is  then  requested  after  a  user 
time  interval  computed  as  follows:  one  whole  second  in  the 
receiver  clock  minus  the  last  computed  offset.  In  other  words, 
after  10,230,000  cycles  minus  the  number  of  cycles 
corresponding  to  the  computed  offset  at  the  previous  epoch. 
The  clock  cycle  is  precisely  the  period  of  the  local 
oscillator,  or  1/10,230,000  =  97.7  ns. 

In  this  way,  the  navigator  is  guaranteed  to  get  the 
next  set  of  measurements  at  the  clock  cycle  transition  that  is 
expected  (from  the  measurements  at  the  previous  epoch)  to  be 
the  closest  to  the  next  even  GPS  second.  When  the  solution 
corresponding  to  the  next  epoch  is  computed,  the  receiver  will 
report,  in  addition  to  the  solution  itself,  how  many  clock 
cycles  the  measurement  epoch  was  shifted,  relative  to  the 
receiver  clock  cycles  counter. 

F.  SOFTWARE  AND  FIRMWARE 

As  a  brand  new  piece  of  hardware,  the  receiver,  as  well  as 
its  supporting  software,  is  subject  to  some  upgrading.  Most  of 
the  research  was  done  on  data  from  receivers  with  a  D020 
firmware  version.  Some  of  them  were  taken  with  an  older 
version,  DOll.  The  major  difference  between  them  is  the  time 
recovery  option,  not  available  in  the  DOll  version.  Data 
formats  are  slightly  different,  too. 

Two  different  CDU  software  versions  were  used,  too: 
version  l.Ol  and  version  1.06.  The  main  difference  is  the 
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ability  of  version  1.06  to  control  the  data  flow  in  the  raw 
data  port.  Unless  specifically  stated,  all  the  information 
here  contained,  refers  to  firmware  version  D020  and  software 
version  1.06. 
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zy.  THE  RAW  DATA  MESSAGE  AMD  THE  COMPUTATION  OF  THE 
MEA8URABLES 


This  chapter  is  devoted  to  the  contents  of  the  raw 
measurement  records,  and  how  all  the  numbers  should  be 
combined  to  build  the  actual  measurements.  The  goal  of  the 
research  is  to  characterize  the  physical  clock  of  the 
receiver.  Because  of  the  high  rates  of  drift  in  the  clock,  the 
receiver  has  been  provided  with  the  Time  Recovery  option.  It 
means  that  the  receiver  can  slew  the  time  of  the  measurement 
to  keep  the  offset,  relative  to  the  GPS  time,  at  the  30  meters 
(100  ns)  level.  That  makes  it  possible  to  take  simultaneous 
measurements  with  different  receivers.  But  if  time  tags  are 
being  adjusted,  all  the  measurements  taken  at  a  given  epoch 
will  contain  the  jump  corresponding  to  that  adjustment.  Let  us 
see  first  how  the  variables  are  reported,  to  see  then  how  the 
measurements  should  be  computed.  Finally,  the  effect  of  the 
time  tag  adjustment  on  the  reported  measurements  will  be 
discussed. 

A.  THE  MX4200  RAW  DATA  MESSAGE 

All  the  quantities  involved  in  measurements,  solution  and 
time  recovery  are  output  by  the  receiver  through  its  raw 
measurement  port.  They  are  reported  in  record  types  1  or  2,  8 
and  969  respectively.  [Ref.  8]  gives  a  description  of  them. 
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1.  Raw  Beasuramant  data  raeord,  typas  i  and  2 

The  raw  measurements  come  out  in  one  of  two  operator's 
selected  format.  Block  I's  have  the  data  for  one  satellite  at 
a  time.  Block  2's  are  compressed  version  of  block  I's  and 
contain  the  information  for  all  the  channels  at  a  time.  Block 
2's  take  about  one  sixth  the  space  as  the  equivalent  data  in 
block  I's.  Although  block  2's  contain  only  printable 
characters,  the  data  is  not  meaningful  till  it  has  gone 
through  an  expansion  algorithm. 

Listed  in  Table  1  is  a  sample  of  a  raw  data  set,  in  a 
block  1  format,  with  all  its  variables  named. 


TABLE  1  RAW  DATA  RECORD  AND  ITS  CONTENT 


ID 

CH 

PRN 

USR  MS 

CHL  MS 

PHI 

CODE  PHFR 

CR 

SNR 

H 

1 

1 

17 

353520500 

353519520 

16478859 

128 

-44 

124 

43 

+ 

1 

2 

20 

353520000 

353519919 

8595657 

114 

-12 

123 

42 

- 

1 

3 

23 

353520000 

353519919 

3547947 

25 

72 

125 

43 

- 

1 

4 

16 

353520000 

353519919 

9408528 

96  -115 

125 

43 

■f 

1 

5 

3 

353520000 

353519918 

20774479 

138  -119 

125 

42 

- 

Name  Deecription 


ID 

CH 

PRN 

USR_MS 

CHL_MS 

PHI 

CODE 

PHFR 

CR 

H 


Record  type,  1 

Channel  number 

PRN  or  eatellite  identifier 

User  time  of  measurement,  milliseconds 

Channel  or  transmit  time,  integer  milliseconds 

Integer  number  of  LI  cycles,  integrated  carrier  phase 

Raw  code  offset,  in  integer  number  of  LI  cycles 

Fractional  part  of  integrated  carrier  phase  in  units  of 

1/256  LI  cycles 

Costas  ratio:  estimate  of  maximum  error  in  phase 
measurement 

Half-cycle  ambiguity  resolution  indicator 
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Several  issues  should  be  pointed  out  here: 

•  Reported  receiving  and  transmit  time  are  in  different  time 
frames:  ^User  time'  and  'GPS  time'  respectively. 

•  The  +/-  sign  at  the  end  of  each  record  ir  described  in 

[Ref  7:p.  4]  as  a  half>cycle  phase  indicator.  This 

resolves  the  half  a  cycle  ambiguity  in  the  phase 
measurement.  When  a  -  sign  is  reported,  half  a  cycle 
should  be  added  to  the  resulting  phase  measurement 

•  For  a  record  type  2,  USR_MS  and  PHI  are  reported  in  a 
slightly  different  way  that  in  block  I's. 

2.  Solution  record/  type  8 

Table  2  shows  a  solution  record  with  all  its  fields 
named  as  the  variables  they  represent. 


TABLE  2.  SOLUTION  RECORD  AND  ITS  CONTENT 


ID  TREC  LAT  LON  HT  VN  VE  VU  VCK  M  N  DP 

8  353513.00  36  35.9357  -121  52.5952  -49.3  0.0  0.0  -0.0  86.6  356 

name  Description 

ID  Record  id . ,  type  8 

TREC  Receiving  time  in  seconds  within  the  week 

LAT  Latitude,  degrees  and  minutes,  WGS  84  Reference  System 

LON  Longitude,  same  as  before 

HT  Height  above  the  ellipsoid,  meters 

VN  North  velocity  component,  meters  per  second 

VE  Bast  velocity  component 

VU  Upwards  velocity  component 

VCK  Clock  drift,  meters  per  second 

M  Navigation  solution  mode  (DR,  2-D,  3-D) 

N  Number  of  satellites  concurrent  in  the  navigation 

solution 

DP  Geometrical  Dilution  of  Precision 


Several  items  to  note  are: 

•  Receiving  time,  TREC,  can  be  optionally  reported  in  UTC  or 
GPS  time.  The  solution  block  following  the  data  that 
generated  it,  will  have  a  time  tag  several  seconds  smaller 
than  the  data  if  UTC  is  selected  (the  default) . 
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•  Solutions  will  be  at  even  seconds  when  time  recovery  (TR) 
is  on.  TREC  will  show  the  clock  offset  if  TR  is  disabled. 
When  a  navigation  solution  cannot  be  computed,  the  TR 
offset  is  held  fixed. 

•  The  number  of  digits  in  each  field  limits  the  resolution 
of  the  navigation  solution:  1/10000  minutes  in  latitude 
and  longitude,  a  decimeter  in  height  and  one  tenth  of  a 
meter  per  second  in  velocity. 

•  CKV  is  the  clock  frequency  drift,  as  computed  by  the 
Kalman  filter,  scaled  to  m/s. 

3.  Time  recovery  record,  type  969 

The  receiver  can  output  a  timing  pulse  every  second, 
when  is  set  to  'time  recovery'.  Under  these  circumstances, 
record  969  provides  useful  information  about  the  estimated 
accuracy  of  that  pulse  and  status  of  the  receiver  clock.  Table 
3  shows  its  contents.  Most  important  for  post  processing,  this 
record  contains  the  information  about  the  adjustment  the 
receiver  has  made  to  its  internal  clock. 


TABLE  3.  TIME  RECOVERY  RECORD  AND  ITS  CONTENT 


ID  CLKMV  FOFF  CLKOFF  ERR 

969  -87.916  -455  -0.000000615  41 


Name  Description 

ID  Record  id.,  type  969 

CLKMV  'Clock  moved'  last  time  of  measurement,  in  meters 

FOFF  Frequency  offset  of  receiver  clock,  relative  to  LI 

carrier  frequency,  cycles 

CLKOFF  Receiver  clock  offset  from  GPS  time 

ERR  Error  in  the  last  time  pulse,  in  units  of  the  local 

oscillator  period  (10.23  Mhz) 


FOFF  (frequency  offset)  corresponds  to  VCK  in  a 
position  record,  expressed  in  cycles  of  the  LI  frequency.  ERR 
(error  in  the  last  time  pulse)  has  to  do  with  hardware  delays 
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in  generating  the  timing  pulse  and  is  not  pertinent  to  this 
analysis. 

CLKOFF  (clock  offset)  is  the  last  estimated  error  in 
the  'corrected'  internal  clock  at  the  last  measurement  epoch. 
When  solutions  are  not  computed,  CLKOFF  stays  at  the  last 
computed  value.  Finally,  CLKMV  ('clock  moves')  is  the  amount 
of  time  the  clock  was  adjusted  prior  to  the  last  measurement. 
It  is  given  in  units  of  length. 

B.  COMPUTATION  OF  THE  MEASURABLES 

Let  us  examine  how  to  compute  the  raw  measurements 
(pseudorange  and  phase-derived  range)  from  the  quantities 
reported  in  the  raw  data  message.  Time  recovery  will  be 
assumed  off.  How  enabling  TR  affects  the  measurement 
computation  will  be  discussed  later. 

The  MX4200  Technical  Manual  provides  an  algorithm  to 
compute  pseudorange  and  phase  range  from  the  reported 
quantities  in  records  type  1  [Ref.  8;p.  C-4].  It  refers  to 
B003  and  D003  firmware  versions: 

=  -{PHI  +  E^)  X  (16) 


^  USRMS 
1000 


(17) 
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CHLMS 

1000 


P*  ^  CODE 


P  =  (T,  -  rj  c 


(18) 

(18) 


Where 

•  X  =  c/f  =  LI  carrier  wavelength 

•  c  =  speed  of  light,  299,792,458  m/s 

•  /  =  LI  carrier  frequency,  1575.42  MHz 

•  T,  =  receiving  user  time 

•  r,  =  transmit  GPS  time 

•  p  =  pseudorange 

•  p^  =  phase-derived  range,  or  biased  phase  range. 

The  April  1991  Firmware  Release  Announcement  [Ref.  9:p.  4] 
adds  a  new  issue  to  the  measurements  computation  for  the 
upgraded  B020,  D020  and  N020  firmware  versions:  half  a  cycle 
should  be  added  to  the  reported  phase,  when  the  half  cycle 
indicator,  H,  is  Then,  reported  phase  measurement  will 

be: 

<p={PHI*  cycles, LI  frequency  (20) 

256  2 

where  a  is  one  when  H  is  0  otherwise. 

Let  t,,  tj,  be  the  reported  USR_MS  and  CHL_MS,  as  expressed 
in  seconds: 
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USRJtS 

1000 


(21) 


CHLMS 

1000 


(22) 


and,  finally,  let  r  be  the  reported  CODE,  in  LI  cycles.  Then, 
expressions  (16)  to  (19)  become: 

P*  =  -0X  (23) 

r,  =  t,  (24) 

T,  *  t,  ^  (25) 

p  =  (T.  -  T,)c  (26) 

Let  us  see  now  that  the  reported  p^  is  in  fact  the  phase 
range,  biased  by  the  integer  milliseconds  t.-t,.  As  will  be 
seen,  the  pseudorange  is  computed  as  the  sum  of  the  phase 
range  and  the  small  value  CODE.  This  is  a  consequence  of  the 
phase  aided  code  tracking  loop,  described  in  Chapter  III. 

The  measured  range  is  in  fact 
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p=At  c 


(27) 


where  At  expresses  the  time  delay  between  transmit  and 
receiving  time  of  a  given  edge  of  the  code  that  modulates  the 
signal.  It  is  measured  in  seconds  of  physical  time.  This  kind 
of  seconds  cannot  be  realized  other  than  by  means  of  a  count 
of  clock  cycles. 

When  the  tracking  is  first  locked,  the  number  of  cycles 
the  replica  was  delayed  can  be  scaled  by  154  (which  is  the 
number  of  LI  cycles  in  a  clock  cycle)  .  So,  the  measured  range 
can  be  expressed  in  terms  of  an  integrated  Doppler  phase, 
i.e.,  carrier  cycles: 

Pq  “  00  ^  (28) 

where  <Po  times  the  number  of  clock  cycles,  and  X  is 
the  carrier  wavelength.  The  error  in  the  first  measurement 
will  be  propagated  to  all  the  measurements  thereafter.  It  is 
the  bias  in  the  phase  derived  range.  The  initial  phase  count 
can  be  split  into  two  parts;  an  integer  number  of 
milliseconds,  by  excess,  minus  a  submillisecond  phase  count: 

~  1000^  ”  ~  (29) 

~  0mj,(A 
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Here  Mg  is  the  initial  integer  nvunber  of  milliseconds,  and  mg 
is  the  same  quantity,  expressed  in  seconds. 

The  output  from  the  correlator  is  used  to  correct  the 
Doppler  frequency  and  time  delay  models,  so  the  process  is 
self  sustained.  The  change  in  range  to  the  satellite  is 
measured  by  integrating  the  Doppler  frequency.  It  is  added  to 
the  previous  <Pg.  Because  the  propagation  meditm  is  dispersive, 
phase  velocity  is  larger  than  group  velocity.  Therefore  a  time 
delay  equivalent  to  that  Doppler  count  will  not  be  enough  to 
get  maximum  output  in  the  correlator.  An  extra  time  delay  is 
introduced,  and  this  is  known,  in  the  Magnavox  notation,  as 
CODE  (F) ,  reported  in  carrier  cycles.  Then,  after  a  certain 
interval,  the  range  can  be  represented  by 


/),.  *  mgC  -  (<t>^  g  +  A0,  r;)X  = 


(30) 


Here  m,  is  an  integer  number  of  milliseconds.  The  range 
measurement  is  thus  divided  into  a  whole  number  of 
milliseconds  and  a  count  of  extra  cycles.  Whenever  0^,  gets  to 
zero  cycles  (or  to  1,575,420  if  the  range  is  increasing),  it 
is  reset  to  1,575,420  (or  to  0).  At  the  same  time,  a  new 
millisecond  is  added  (or  subtracted)  to  m,.  Figure  3  shows  how 
those  milliseconds  adjustments  in  phase  are  done.  The  jagged 
curve  at  the  bottom  is  what  the  receiver  reports  as  phase, 
scaled  to  range.  If  that  reported  phase  range  is  subtracted 
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Figure  3.  Millisecond  adjustments  in  the  reported  phase 

from  the  range  corresponding  to  the  integer  number  of 
milliseconds,  the  uppermost  curve,  phase-derived  range,  is 
obtained.  Thus,  eg.  (16)  and  (23)  should  be  properly  adjusted 
to  have  the  actual  phase  derived  range: 

P*  =  (t.  -  t,)c  -  (31) 


If  in  addition,  the  CODE  range  equivalent  is  subtracted 
from  the  phase  derived  range,  the  pseudorange  is  obtained, 
i.e.,  the  p  curve  in  Figure  3. 
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Code 


r  =  -  p)  (32) 


gives  the  difference  between  measured  phase  range  and 
pseudorange.  The  hardware  architecture,  thus,  allows  the  code 
component  of  the  loop  to  track  only  changes  in  the  difference 
between  biased  phase  range  and  pseudorange.  This  difference 
consists  of: 

•  Ionospheric  delay,  that  enters  with  opposite  sign  in  the 
phase  measurement  and  .ange  (time  delay)  measurement. 

•  Multipath  effects,  that  are  different  for  phase  and  code, 
too.  Their  interference  patterns  are  not  alike  because  the 
difference  in  wavelengths  of  carrier  (0.19  m)  and  C/A  code 
(300  m) . 

•  Difference  in  the  noise  level  of  pseudorange  and  phase 
measurements.  It  is  about  two  orders  of  magnitude  larger 
for  pseudorange. 

•  Finally,  it  contains  the  original  bias  in  the  phase  range. 

These  effects,  combined  all  together,  have  a  much 
smaller  dynamics  than  the  actual  range  to  the  SV.  Thus,  it  is 
possible  to  keep  it  tracked  in  a  much  narrower  bandwidth. 

C.  RAW  VERSUS  REPORTED  NEASUREMEMTS 

Let  us  see  now  how  to  set  the  receiver  to  time  recovery 
option  affects  the  reported  measurements,  so  we  can  adjust 
them. 


48 


1.  Thr««  different  tine  f rases 


In  the  literature  dealing  with  precise  time,  two 
different  concepts  are  common:  'physical'  clocks  as  opposed  to 
'paper'  clocks.  Let  as  assume  we  have  a  wall  clock  to  watch, 
that  lags  an  indeterminate  number  of  seconds  in  a  certain 
period.  Current  time  is  going  to  be  requested  once  in  a  while. 
To  deal  with  that  lag,  we  carefully  listen  to  the  broadcast 
timing  signals  in  WWV.  At  the  same  time,  we  keep  a  record  of 
the  number  of  even  seconds  our  wall  clock  lags.  Vfhenever  time 
is  requested,  we  will  report  what  our  wall,  physical  clock 
says,  plus  the  number  of  seconds  we  know  for  sure  it  lagged 
last  time  the  record  was  updated.  This  whole  process  is  known 
as  a  'paper  clock'.  Neither  of  those  two  times  is  true  time, 
which  is  only  known  at  timing  signal  epochs.  Then,  three 
different  time  frames  are  present  in  the  example  above. 

There  are  three  different  time  frames  in  the  receiver. 


too; 

•  Absolute  time,  corresponding  to  GPS  time  in  this  case.  It 
is  only  estimated  whenever  the  receiver  gets  a  navigation 
solution.  The  estimate  has  enough  accuracy  to  decide  how 
many  cycles  the  measurement  epoch  lags  the  GPS  time.  This 
corresponds  to  the  absolute  time  in  the  example  above. 
Absolute  time  will  be  represented  by  T^. 

•  Physical  time,  that  is  kept  by  the  receiver  physical 
clock;  it  is  driven  by  a  freewheeling  oscillator,  at  a 
nominal  frequency  10.23  Mhz.  It  is  represented  by  T^. 

•  Paper  time:  this  is  the  time  the  receiver  reports  in  the 
data  stream,  as  the  time  of  the  measurement.  This  is  known 
as  user  time  and  will  be  expressed  as  T^. 
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Measurement  epochs  are  at  even  seconds  in  user  time. 
But  user  time  is  an  artifact.  The  receiver  needs  to  compute 
when,  in  its  physical  clock  time,  the  measurements  are  to  be 
taken.  That  is,  precisely  at  which  transition  of  the  physical 
clock  counter. 

The  cycle  of  measurement  and  time  adjustment  begins 
with  a  navigation  solution.  The  estimate  of  the  clock  offset 
and  drift  (from  navigation  solution)  are  used  to  predict  the 
number  of  physical  clock  counts  between  the  last  measurement 
and  the  next  true  GPS  second  mark.  The  difference  between  this 
value  and  the  nominal  count  of  Nq  =  10,230,000  (one  second  at 
10.23  MHz),  is  the  value  S  reported  in  the  969  record,  i.e., 
CLKMV  {S  in  fact  has  been  scaled  to  range,  to  be  reported  as 
CLKMV)  .  When  the  last  measurement  was  taken,  the  cycles  count 
was  reset  to  zero.  The  next  measurement  will  be  taken  when  the 
cycle  counter  of  its  physical  clock  gets  to  +  S,  i.e.,  at 
a  physical  time 


m  =7* 

11,1  *1  u,i 


No  + 


(33) 


The  same  process  occurred  in  the  previous  measurement  epoch, 
so  this  number  of  cycles  has  been  added  to  the  previous  clock 
shift.  This  is 
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(34) 


and 


=  T,  *  a  *  1)  *  ‘ 

^  A 


(35) 


Because  no  offset  existed  at  initial  tine 


=  *  i 


(36) 


SO 


~  ^p4*l 


(37) 


Figure  4  shows  the  relationship  between  these  three 
time  franes.  T,  is  the  true  GPS  time;  the  interval  between  T^, 
and  is  precisely  10,230,000  times  the  nominal  period  of 
the  oscillator.  Once  the  receiver  has  a  measurement,  at  T,, 
(which  is  an  even  second  in  user  time) ,  the  time  and  frequency 
offset  in  the  navigation  solution  are  used  to  compute  6,,  i.e. , 
the  number  of  cycles  for  the  next  measurement  epoch  to  be 
shifted. 

A  new  cycles  count  starts  then  and  the  next 
measurement  will  be  taken  after  Ng  +  cycles.  Cycles  in  the 
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ith  GPS  second  GPS  second 


Figure  4.  Time  frames  in  the  receiver.  T^,  Ty  and  Tp  are 
'GPS',  'user',  and  'physical'  time  respectively.  5;  is  the 
clock  shift  at  epoch  i+1.  Aj  is  the  integrated  clock  shift. 


physical  clock  do  not  have  precisely  the  same  period,  though. 
When  that  count  gets  to  the  end,  and  a  new  measurement  is 
taken,  a  different  estimate  is  found.  The  next  measurement 
will  be  taken  after  cycles  of  the  physical  clock. 

2.  Raw  Baaauremeats  computation 

The  noise  in  the  physical  clock  of  the  receiver  is  the 
subject  of  the  research.  In  a  Time  Recovery  off  situation  no 
clock  shifting  is  done.  Thus,  both  measurements  and  time 
tagging  are  expressed  in  terms  of  the  physical  clock  readings. 
User  and  physical  times  are  synonymous. 

Contrarily,  if  Time  Recovery  is  on,  the  time  of  the 
measurement  i  has  been  shifted  by  A,.;  cycles,  with  respect  to 
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the  physical  time.  The  time  tagging  process  is  driven  by  the 
physical  clock.  The  timing  error  is  the  error  in  the  physical 
clock.  Thus,  to  compute  it,  we  need  the  reading  of  the 
physical  clock  to  which  A,.,  has  been  added  to  get  the 
measurement  epoch,  i.e. 


T 

p.> 


(38) 


Now,  if  the  measurement  epoch  has  been  delayed  until 
T^j,  the  measurements  at  that  epoch,  as  expressed  in  terms  of 
user  time,  have  been  biased  by  the  range  corresponding  to 
cycles: 


-  T,)c  -  ip\ 

=  ^  -  T,)c  -  0X 

•‘o 


(39) 


Thus,  to  retrieve  the  raw  measurements,  the  same 
adjustment  should  be  done: 


(T.  -  r,)c  -  0X 


(40) 


and  the  same  applies  to  the  pseudorange.  By  comparing  these 
ranges  with  the  true  ranges  it  is  possible  to  see  the  error  in 
the  physical  clock. 
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Now,  raw  measurements  can  be  obtained  in  terms  of  the 
reported  quantities  as  follows: 


/>4»( 


1000 


(41) 


p=p^-CODE  X 


regardless  the  Time  Recovery  option  that  was  selected.  The 
physical  time  corresponding  to  that  epoch  will  be: 


_  USRJiS 
”  1000 


i-i 

-  J^CLKMVj 


(43) 
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V.  CHARACTERISTICS  OP  THE  RECEIVER  CLOCK 


From  the  information  contained  in  Chapter  IV,  it  is 
possible  to  retrieve  the  raw  measurements  the  receiver 
performs.  They  are  contaminated  by  many  different  biases  and 
errors.  Those  external  to  the  receiver  are  not  the  main  issue 
here.  They  should  be  removed  from  the  measurements,  prior  to 
the  analysis,  though.  It  is  mainly  the  clock  characteristics 
which  determines  the  accuracy  in  the  measurement.  This  will  be 
the  subject  of  the  present  and  next  chapters. 

A.  DETERMINISTIC  EFFECTS  PRESENT  IN  THE  MEASUREMENT. 

All  those  effects  present  in  the  measurement  that  can  be 
somehow  modeled,  will  be  considered  as  deterministic. 
Stochastic  noise,  on  the  other  hand,  cannot  be  modeled.  Their 
sources  may  be  identifiable  as  either  noise  or  unremoved 
biases.  They  can  only  be  quantified  as  a  bulk. 

1.  External  to  the  clock 

All  those  biases,  mentioned  in  Chapter  II,  as  due  to 
the  satellite  (clock  bias,  position  and  velocity)  and  the 
propagation  medium  (ionospheric  and  tropospheric  delays),  are 
considered  here  as  deterministic  effects  external  to  the 
receiver  clock. 
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2.  Internal  to  tha  clock 


The  accuracy  of  a  clock  is  the  deviation  of  its 
reading  from  an  absolute  time.  In  a  general  sense,  this 
accuracy  can  be  expressed  as  a  polynomial  of  a  certain  degree 
in  At,  plus  a  stochastic  noise  component. 

That  polynomial,  that  models  the  bulk  of  the  clock 
accuracy,  is  considered  as  the  deterministic  part  of  the  clock 
error.  The  order  of  the  polynomial  needed  is  highly  dependent 
on  the  stability  of  the  clock  and  the  elapsed  time  since 
synchronization.  As  a  matter  of  experience,  up  to  a  third 
degree  polynomial  model  should  be  attempted  for  this 
particular  clock,  and  half  a  day  span  of  time.  A  guadratic 
model  would  be  enough  for  periods  shorter  than  one  hour.  Later 
on  it  will  be  seen  that  some  effects  of  temperature  change  can 
be  present.  Its  effect  is  very  hard  to  separate  from  those  due 
to  the  physical  characteristics  of  the  particular  crystal  in 
question. 

In  addition  to  that,  Hellwig  [Ref.  10:p.  27-30]  gives 
four  main  deterministic  effects  in  crystal  oscillators.  Only 
the  two  first  were  studied  here: 

•  Temperature. 

•  Transient  regimes. 

•  Gravitational  forces  and  accelerations. 

•  Line  voltage  variations. 


56 


B.  THE  CLOCK  MODEL 

As  a  fundamental  analysis  tool,  a  FORTRAN  program  was 
written  to  calculate  a  model  of  the  receiver  clock  offset, 
relative  to  GPS  time.  It  is  based  on  the  recorded  raw 
measurements,  in  a  given  observation  session.  Appendix  A, 
contains  a  description  of  it. 

This  program  calculates  a  model  of  the  deterministic  part 
of  the  clock  behavior.  In  addition,  a  further  characterization 
of  the  stochastic  noise  is  possible.  It  is  done  by  analyzing 
the  residuals  of  the  raw  measurements,  from  that  model.  This 
will  be  the  subject  of  Chapter  Vi.  Underlying  this  procedural 
scheme  there  is  an  assumption:  all  the  deterministic  effects, 
other  than  the  clock  stability,  are  not  present  in  the 
measurements.  To  achieve  this,  three  different  approaches  can 
be  taken,  according  to  each  particular  effect: 

•  To  fix  the  conditions  under  which  the  measurements  were 
taken,  so  some  of  the  effects  could  be  considered  as 
negligible.  In  particular:  receivers  at  rest,  fed  by  a 
well  stabilized  power  supply,  well  off  transient  regimes, 
and  carefully  watched  temperature  variations. 

•  To  subtract  from  the  measurements  a  model  of  every 
identifiable  (and  modellable)  effect:  true  range  to  the 
corresponding  SV  and  a  Hopfield  tropospheric  delay  model. 

•  For  ionospheric  delay,  whose  models  are  not  very  reliable, 
just  assume  its  effect  as  negligible.  Spread  over  the 
entire  visible  constellation  of  SV's,  it  introduces  a 
small,  unknown  bias  into  the  clock  offset.  No  significant 
effect  is  introduced  in  the  frequency  offset;  it  can  be 
considered  as  a  secular  variation  over  a  SV  pass. 

Not  meeting  any  of  these  conditions  may  lead  to  a  model  or 

noise  characteristics  dominated  by  the  particular  condition 
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lacking.  On  the  other  hand,  modifying  only  one  of  these 
conditions  allows  its  effect  on  the  clock  to  be  estimated. 
Experiments  on  temperature  variations,  the  effect  of  ephemeris 
changes  and  SA  activation  on  the  noise,  and  duration  of  the 
transient  regime  were  done.  They  will  be  discussed  later. 

A  description  of  the  program  MXCKMO  is  given  here  to  show 
how  all  the  issues  above  were  taken  into  account.  It  is  linked 
to  several  routines  in  the  GPS  packages  FIC5,  GAPP5  and  GENS 
by  the  Applied  Research  Laboratory,  University  of  Texas  at 
Austin.  They  are  listed  in  Appendix  B,  and  a  brief  description 
of  them  will  be  given.  Several  other  low  level  routines,  such 
as  data  decompression,  extraction  of  ephemeris  and  data 
logging,  were  written  by  Dr.  James  R.  Clynch.  They  are  used 
here  as  interfaces  between  the  receiver  output  and  the  main 
analysis  program.  They  are  listed  in  Appendix  C. 

1 .  Data  logging 

An  initial  problem  was  found  in  the  raw  data  stream. 
Clock  shifts  are  given  in  instantaneous  values,  at 
measurements  epochs.  As  shown  in  eg.  (41),  integrated  values 
are  required  to  retrieve  the  raw  measurements. 

To  avoid  logging  all  data,  for  sessions  lasting  a  few 
hours,  a  new  routine  was  needed.  MXR  adds  all  the  reported 
clock  shifts,  up  to  the  current  epoch.  The  result  is  saved,  as 
integrated  clock  shift,  in  a  new  field  added  to  the  969 
record.  In  addition,  MXR  selects  entire  blocks  of  1/2,  8  and 
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969  records  to  be  logged,  at  a  user's  choice  epoch  interval. 
Other  useful  records,  such  as  ephemeris  and  almanac,  can  be 
recorded,  too,  taking  advantage  of  the  user's  ability  to 
request  their  output  from  the  CDU. 

2 .  Data  interface 

Once  the  raw  data  are  recorded,  the  MX42RTN  package 
performs  several  operations  on  them.  One  routine  allows  a 
sequential  extraction  in  entire  blocks  of  raw  data,  navigation 
solution  and  time  recovery,  corresponding  to  the  same  epoch. 
Among  those  operations,  decompression  of  records  type  2  and 
ephemeris  and  almanac  records  are  key  items.  Ephem^ris  data 
can  be  converted  to  FIC  binary  format,  and  saved  in  a  file  by 
means  of  the  program  MX2EPH.  Several  routines  in  the  FIC5  and 
GAPP5  packages  (LOADTAB,  GETEPH,  FINDEPH) ,  handle  them  in  an 
efficient  way.  The  whole  ephemeris  set  is  loaded  in  a  20x3x130 
array.  The  desired  items  are  fetched  in  RAM,  rather  than  in 
files. 

3 .  Input  interface 

This  allows  the  user  to  enter  all  the  parameters  for 
the  computation.  In  particular,  the  following  are  prompted: 

•  Names  of  data,  ephemeris  and  output  files. 

•  Initial  and  final  epochs  for  the  computation  to  be 
completed. 

•  Reference  time  of  epoch  and  degree  of  the  polynomial 
model. 
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•  Specification  of  PRN's  to  be  removed  from  the  model 
computation.  It  can  be  useful  whenever  a  given  SV  inputs 
undesired  effects  in  the  computation. 

•  Specification  of  PRN's  for  which  output  files  are  not 
desired.  MS  FORTRAN  5  imposes  a  limitation  in  the  number 
of  files  can  held  simultaneously  open. 

•  Selection  of  the  ground  position  from  those  stored  in  the 
STATIONS. FIL  file.  LOAOSTA  routine,  in  6APP5  package, 
allows  the  retrieval,  coordinates  transformation  and 
storage  in  memory,  of  the  antenna  position. 

•  Surface  meteorological  data,  to  be  used  in  the  ANLTROP 
function  (GAPP5) .  This  is  a  Hopfield  tropospheric  model. 

4 .  Main  algorithm 

The  algorithm  is  based  on  the  computation  of  the 
coefficients  of  a  polynomial.  It  models  the  physical  clock 
offset,  relative  to  GPS  time.  At  a  given  epoch  i,  that  offset 
is  first  estimated  by  comparing  measured  pseudorange  to  the 
satellite,  with  the  true  range.  The  same  process  is  done  for 
all  the  available  satellites  at  that  epoch  (except  those 
specifically  excluded  from  the  model  computation) .  Between 
initial  and  final  epochs,  all  the  range  residuals  are 
computed.  These  observed-minus-computed  ranges  (OMC's),  scaled 
to  time,  enter  in  a  batch  Ordinary  Mean  Least  Squares 
Adjustment  (OMLS)  to  find  the  coefficients  of  the  polynomial 
that  better  fits  them.  The  adjustement  process  is  described  in 
Appendix  D. 

Five  main  components  can  be  identified  within  the  main 
algorithm:  Raw  data  input,  receiver  clock  correction,  OMC 
computation,  computation  of  the  deterministic  model  and  statistics. 
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a.  Raw  data  input 

These  data  are  read  by  calling  MXRDRW  routine.  In 
turn,  it  loads  an  entire  block  of  data  by  means  of  NXMX42  and 
M42SET  (MX42RTN  package) .  Raw  measurements,  uncorrected  by 
clock  shift,  and  the  clock  shift  itself,  are  returned. 

For  data  logged  in  TR  option,  the  integrated  clock 
shift  at  the  previous  epoch  is  needed.  The  routine  selects 
integrated  clock  shift  from  either  the  previous  or  the  current 
block,  depending  on  the  format  of  the  measurements  record.  If 
raw  measurements  are  in  compressed  format,  the  integrated 
clock  shift  in  the  previous  block  is  needed.  It  is  computed  by 
subtracting  the  instantaneous  clock  shift  from  the  integrated 
clock  shift  in  the  current  block.  If  raw  measurements  are  in 
ASCII  format,  the  integrated  clock  shift  at  the  previous  epoch 
is  in  the  current  block. 

Once  the  main  program  has  the  unadjusted  (by  clock 
shift)  ranges  and  the  clock  shift,  the  raw  ranges  are  computed 
following  eg.  (41)  and  (42) .  The  reading  of  its  physical 
clock,  at  that  epoch,  is  also  computed; 


If  TR  was  off,  or  navigation  solution  failed, 
is  0 .  The  measurements  are  taken  at  an  even  second  of  the 
physical  clock  for  that  case. 
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b.  RBCBiver  clock  correction 

The  physical  time  of  the  receiver  clock  should  be 
now  corrected  by  its  offset,  relative  to  GPS  time.  This  is 
done  by  means  of  the  clock  offset  model  computed  in  the 
previous  iteration.  The  model 

T  (t)  =a+b(t-to)  +1  (t-to)2+|  (t-to)’  (45) 

is  expressed  in  terms  of  a  GPS  time  interval.  It  might  be 
known  only  to  within  a  few  tenths  of  a  second.  But  t  (GPS 
time)  is 

t=t^-T(t)  (46) 

Then,  t  is  first  computed  as 

T‘=a+b(Vto)-^|(Vto)2-H|(Vto)’  (47) 

and  then  iterated  as 

T-=a+b(t^-T-'-to)+£(VT'->-to)'+|(t^-T->-to)’  (48) 

until 

|T‘-T-'i<10-'^ 
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Then  pseudorange  and  phase  range,  as  well  as 


physical  time,  are  corrected  by  this  clock  offset: 

(50) 


P*c=P^-‘rc 


(51) 


(52) 


c.  OMC  computation 

To  perform  the  computation  of  the  true  range,  the 
closest  ephemeris  set  to  the  current  epoch  is  chosen.  This  is 
not  a  requirement  for  data  sets  spanning  less  than  three 
hours.  The  elements  of  the  broadcast  orbit  are  claimed  to  be 
good  within  six  meters  and  valid  for  four  hours.  They  are 
updated  every  two  hours. 

For  periods  longer  than  three  hours,  the  SV 
position  must  be  computed  with  different  sets  of  ephemeris.  It 
has  the  effect  of  introducing  spurious  noise  when  switching 
from  one  to  another.  The  difference  between  ranges,  as 
computed  with  two  consecutive  ephemeris  sets,  is  taken.  This 
difference  is  then  added  to  the  new  range,  so  the 
discontinuity  is  saved. 

To  compute  the  true  range  to  the  SV,  the  GPS  time 
of  the  measurement  is  needed.  Again,  this  is  only  known 
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approximately.  Residuals  off  the  clock  model  have  been 
observed  to  be  as  much  as  100  Km  ,  i.e.,  0.3  ms.  The  SV  moves 
about  4  Km/s.  At  that  velocity,  a  0.3  ms  error  in  time  makes 
an  error  in  the  computed  position  of  about  one  meter. 

Actual  GPS  receiving  time  is: 

t=V(T+€)=t^-e  (53) 


where  e  is  the  residual  from  the  clock  offset  model  r.  The 
range  to  the  SV  is  first  computed  at 

t‘=t,,  (54) 

Tropospheric  model  correction  is  then  applied  to 
this  computed  range,  so 

(55) 


Then,  a  first  estimate  of  the  OMC  and  the  time  residual  is 
computed : 


I  Pc~^c 
€  - - — 


(56) 


This  first  residual  is  then  added  to  the  clock 
offset  model.  A  better  estimate  of  the  true  GPS  time  is  then 
obtained. 
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t2=t‘-€'=tp-(T+e‘)  (57) 

and  new  true  range  and  residual  are  computed.  The  process  is 
iterated  until  the  change  in  the  residual  is  less  than  10'*^ 
seconds.  The  OMC  so  computed 


c 

is  in  fact  the  time  residual  corresponding  to  the  given  SV,  at 
current  epoch  i,  for  the  present  iteration  of  the  adjustment 
process . 

The  routine  to  compute  the  theoretical  range 
follows  an  iteration  scheme,  too.  It  first  computes  the  SV 
position  and  range  for  the  input  time,  i.e.,  receiving  time. 
Then,  it  computes  transmit  time  by  subtracting  the  time  of 
flight,  and  corrects  the  ground  position  by  Earth  rotation. 
The  process  is  iterated  until  convergence.  Relativistic 
effects  are  corrected  for,  too. 

d.  Computation  of  the  deterministic  model 

As  residuals  for  each  SV  and  epoch  are  being 
computed,  a  new  observation  equation  is  added  to  the  system. 
The  IfH  and  //*€  are  being  updated  with  new  entries  every  time 
a  new  measurement  is  processed.  For  N  observation  equations, 
the  size  of  the  elements  in  IfH  ranges  from  N  to  NxAt^.  For 
three  hours  worth  of  data,  at  one  second  logging  rate,  those 
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numbers  are  in  the  order  of  10^  and  10**.  ItH  is  ill  conditioned 
for  inversion  purposes. 

The  normal  equations  system 

7  C59) 


is  solved  by  LU  decomposition  [Ref.  ll:p.  31-43].  Routines 
LUDCMD  and  LUBKSD,  modified  to  double  precision,  were  used. 

The  first  solution  to  the  system,  is  iterated  as 


H‘H  x-H‘H  Xq=WH  dx‘= 
/f'  €-H‘H  Xo 


(60) 


SO  the  system  can  be  solved  for  Ar  now,  and 

x=Xo+A?f  (61) 

Iterations  are  done  until  the  change  in  the 
parameters  is  less  than  where  j  is  the  current 

iteration  of  the  adjustment  process.  This  guarantees  a  change 
less  than  10'^  in  the  fourth  iteration  of  the  main  process.  For 
three  hours  worth  of  data,  at  1  second  rate,  the  change  in  the 
cubic  term,  worst  case,  would  be  less  than  (5  x  10^)*xl0‘^,  in 
the  order  of  10'**  seconds. 

The  convergence  criteria  for  the  adjustment  process 
was  set  to  halt  the  process  when  the  variation  in  any 
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parameter  would  produce  a  clock  offset  increment  lesser  than 
1  ns,  worst  case. 

e.  Statistics  and  edition  of  data 

The  raw  data  file  might  contain  several  anomalous 
records.  They  degrade  the  model  computation.  The  main  routine 
implements  different  algorithms  to  edit  them  out  of  the 
computation.  Each  computed  OMC  is  tested  against  a  window.  The 
size  of  that  window  is  reset  to  five  times  the  a-posteriori 
variance  of  the  residuals  in  the  previous  iteration. 
Pseudoranges  failing  this  test  are  excluded  from  the  model 
computation. 

Statistics  on  the  residuals  after  adjustment  are 
given  on  a  SV  basis.  It  allows  the  user  to  exclude  those  SV's 
that  show  anomalous  behavior,  from  a  later  model  computation. 
Statistics  on  Doppler  measurements  also  are  done  on  a  SV 
basis.  Doppler  measurements  do  not  enter  the  model 
computation.  Increment  in  phase  range,  between  previous  and 
current  epoch,  is  computed  every  epoch.  This  difference, 
divided  by  the  computed  interval,  gives  the  measured  Doppler. 
The  increment  in  true  range,  divided  by  the  same  interval, 
gives  the  true  Doppler  shift,  corresponding  to  the  true  change 
of  the  SV  position.  The  difference  between  them  is  called 
Doppler  residual.  It  is,  in  fact,  the  frequency  offset  in  the 
receiver  clock,  at  the  measurement  epoch. 
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Some  of  the  phase  range  measurements  can  be  wrong. 
Then,  both  the  current  and  next  Doppler  residuals  will  be  in 
error.  To  avoid  that  problem,  a  dynamic  window  edits  those 
anomalous  measurements  out  of  the  statistics  computation.  The 
size  of  the  window  is  reset  once  a  new  Doppler  residual  is 
computed.  It  is  designed  as  a  recursive  filter  of  equation: 

y,=0.1  x,.  +  0.9  y,.,  (62) 

where  y,,  x,  are  the  smoothed  and  actual  residuals.  A  rejected 
Doppler  residual  is  flagged  in  the  output  file. 

The  way  that  record  is  flagged  depends  on  the  kind 
of  anomaly  found:  acquisition,  reacquisition,  gap  in  the  data 
or  anomalous  phase  range.  At  a  later  stage  in  the  analysis, 
that  flag  gives  an  indication  on  the  confidence  that  should  be 
credited  to  that  measurement.  This  is  helpful  to  decide 
whether  or  not  to  use  that  particular  record  for  noise 
analysis  purposes. 

5 .  Outputs 

One  residuals  file  per  selected  SV  is  output.  It 
contains  the  following  items: 

•  Modeled  GPS  time  of  the  measurement  epoch. 

•  Phase-derived  range  residuals. 

•  Doppler  residuals. 

•  Doppler  residuals  flag. 
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READ  10276  DATA  RECORDS  IN  1  TH  LOOP 

READ  10276  DATA  RECORDS  IN  2  TH  LOOP 

READ  10276  DATA  RECORDS  IN  3  TH  LOOP 

prn  12  switched  to  TOE  345600.  eph.  itm.  at  336900. 

prn  13  switched  to  TOE  345600.  eph.  itm.  at  336900. 

prn  15  switched  to  TOE  345600.  eph.  itm.  at  336900. 

prn  2  switched  to  TOE  345600.  eph.  itm.  at  336900. 

prn  6  switched  to  TOE  345600.  eph.  itm.  at  336900. 

READ  10276  DATA  RECORDS  IN  4  TH  LOOP 


PSEUDORANGE 

RESIDUAL 

STATISTICS 

PER  SATELLITE  : 

PRN 

N. 

POINTS 

AVG. 

STD  DEV 

12 

2042 

5.81 

16.672 

13 

2096 

-5.79 

15.611 

15 

1946 

-.09 

15.974 

2 

2096 

-3.73 

15.047 

6 

2096 

3.96 

15.849 

COMPUTATION  COMPLETED  BETWEEN  336900  AND  339000  USING  5  SATELLITES: 
12,  13,  15,  2,  6, 

4  ITERATIONS  WERE  DONE  ON  DATA  FILE  m22115.149 

3  PARAMETERS  CLOCK  MODEL  AT  TIME  OF  EPOCH  WN  594.  SOW  337950.: 
-.144134995E-02  58027»412E-06  . 756479538E-12  . OOOOOOOOOE+00 

A  POSTERIORI  VARIANCE  OF  UNIT  WEIGHT:  .300662E-14  DET:  .108760E+29 

RMS  POST  FIT  RESIDUAL  .055  MICROSECONDS,  16.44  METERS 

DOPPLER  RESIDUALS  USING  10271  POINTS:  AV6  .03020  ,  SD  .0815  M/S 


rigurc  5.  A  MXCKMO  summary  file  printout. 


A  summary  file  is  output,  too.  A  sample  printout  is 
shown  in  Figure  5.  It  contains  the  following  items: 

•  Coefficients  of  the  polynomial  clock  model  and  its 
reference  time. 

•  Mean  and  standard  deviation  of  pseudorange  and  Doppler 
residuals  as  a  bulk. 

•  Mean  and  standard  deviation  of  pseudorange  residuals, 
itemized  per  SV.  A  certain  SV  can  thus  be  rejected  from  a 
later  computation,  based  on  these  estimators. 

•  Epoch  when  a  new  ephemeris  set  was  used  for  the  first 
time,  for  a  given  SV,  as  explained  in  4.c.  This  allows  the 
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user  to  correlate  certain  noises  with  those  changes  in 
ephemeris. 

C.  LONG  TERM  DETERMINISTIC  EFFECTS 

Two  major  deterministic  effects,  among  those  mentioned  in 
V.A.2,  have  been  analyzed:  temperature  changes  and  transient 
regimes.  Both  of  them  have  the  same  origin:  the  variations  of 
the  temperature  in  the  oscillator.  But  they  affect  the 
operational  characteristics  of  the  receiver  in  different  ways. 
Both  of  them  affect  the  frequency  stability  of  the  clock,  in 
the  range  of  hours.  To  correlate  frequency  variations  with 
changes  in  temperature  and  warm-up  times  shows  the  effect. 
Unfortunately,  the  experiments  were  done  under  conditions  that 
can  hardly  be  considered  as  laboratory- like.  This  is  not  a 
major  limitation,  though.  The  MX4200  is  not  a  high-performance 
device.  Each  particular  oscillator  is  likely  to  have  different 
characteristics.  At  least,  this  is  the  case  for  the  two 
receivers  tested  here.  Some  trends  and  ranges  can  be  inferred 
from  the  results  of  these  tests,  although  generalization  of 
the  results  to  all  the  receivers  belonging  to  the  same  family 
is  not  warranted  from  this  sample  of  two. 

1.  Transient  regimes 

An  observation  session  was  held  on  June  18th,  1991 
(day  169)  under  the  following  circumstances: 

•  Receiver  model  707361-803,  SN  2611,  and  antenna  model 
819760-801,  SN  0391.  Hereafter,  this  set  will  be 
considered  as  'receiver  2'. 
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•  TR  on,  static  mode. 

•  The  receiver  was  off  at  least  24  hours  before  the  session. 

•  Data  were  logged  between  16.00  PST  (23.00  GMT)  and  19.30 
PST  (02.30,  June  19th,  GMT)  at  one  second  re"e. 

•  Between  3  and  5  SV's  were  tracked  along  the  session. 

•  The  receiver  was  inside  a  small  building  whose  door  was 
opened  and  closed  several  times  during  the  session.  Short 
term  temperature  changes,  up  to  two  degrees,  can  be 
estimated . 

A  linear  model  of  the  clock  offset  was  generated, 
using  MXCKMO.  Doppler  residuals,  off  that  model,  were 
extracted.  The  uppermost  part  of  Figure  6  shows  the  residuals 
plotted  against  time.  Residuals  off  a  linear  model  allows 
estimating,  at  a  glance,  the  length  of  the  transient  regime. 
Doppler  residuals  in  fact  represent  the  frequency  deviation 
off  the  modeled  frequency  offset. 

A  transient  regime  shows  up  prominently  within  the  two 
first  hours.  The  jagging  in  the  curve  could  be  due  to  short 
term  temperature  variations  in  the  room.  Its  effect  in  the 
frequency  offset  will  be  seen  in  the  following  paragraph. 

As  a  matter  of  generality,  a  second  example  is  shown 
here.  The  observation  session  was  held  on  April  26th,  1991 
(day  116) ,  as  follows: 

•  Receiver  model  707361-803,  SN  1780,  and  antenna  model 
818766-801,  SN  0391.  This  set  will  be  considered  as 
'receiver  1'  hereafter. 

•  Firmware  version  DOll  was  still  installed.  Thus,  TR  option 
was  not  available. 

•  Receiver  was  off  at  least  24  hours  before  the  experiment. 
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•  Data  were  logged  between  15.00  PST  (22.00  GMT)  and  18.30 
PST  (01.30,  April  27th,  GMT)  at  60  seconds  rate. 

•  Four  to  six  satellites  were  tracked. 

•  It  is  presumed  to  have  been  exposed  to  a  fairly  stable 
ambient  temperature. 

The  temperature  variation,  in  this  case,  is  estimated 
to  be  even  less  than  before.  The  receiver  was  inside  a  fairly 
well  insulated  room.  An  external  heat  source  (TI  4100  GPS 
receiver)  was  present  in  the  room.  No  major  air  exchange,  in 
a  short  period,  happened. 

The  same  process  as  above  was  followed,  and  the 
Doppler  residuals  are  shown  at  the  bottom  of  Figure  6.  A 
transient  regime  comparable  to  the  one  above  can  be  seen 
during  the  first  two  hours.  Different  ambient  temperatures 
between  sessions  must  count  for  the  initial  difference  in  the 
size  of  the  residuals.  Different  performances  in  both 
oscillators  should  be  expected,  too. 

2.  Temperature  effects 

This  is  a  major  issue  in  the  performance  of  the 
receiver.  It  was  recognized  in  trying  to  model  a  long  term 
data  set.  Figure  7  shows  the  Doppler  residuals  off  a  linear 
model  generated  using  the  data  between  18  and  34  hours.  It 
shows  a  well  developed  diurnal  pattern.  One  first  thinks  of 
the  temperature  as  its  cause.  To  confirm  it,  a  thermal  shock 
experiment  was  conducted  lasting  six  hours. 
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Olurnal  uaLriation  In  th«  frequency 
Racsiuar  n\-,  residual:  off  a  linear  model 


At  a  first  stage,  receivers  1  and  2  were  set  to  log 
data,  at  ambient  temperature,  for  two  hours.  Transient  regimes 
were  avoided  by  turning  the  receivers  on  days  ahead. 
Temperature  was  recorded  with  a  meteorologic  thermograph.  Its 
readings  were  tested  against  a  calibrated  mercury  thermometer. 
An  agreement  was  found,  within  the  resolution  of  the  graph. 

In  a  second  stage,  both  receivers,  thermograph  and 
mercury  thermometer  were  put  inside  a  refrigerator.  It  was 
turned  on  several  hours  in  advance.  The  lid  was  sealed  with 
foam  to  prevent  heat  exchange.  This  stage  lasted  for  two  more 
hours.  Then,  all  the  equipment  was  removed  from  the 
refrigerator.  The  ambient  temperature  in  the  room  was  close  to 
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30  centigrade.  Two  more  hours  of  data  were  collected,  so  the 
receivers  could  reach  a  steady  state. 

Data  logged  every  six  seconds  were  used  to  get  a 
linear  model  of  the  clock  offset  in  both  receivers.  Doppler 
residuals  off  their  models  are  plotted  in  Figure  8,  along  with 
temperature.  A  certain  thermal  inertia  should  be  expected  in 
the  thermograph.  It  was  in  fact  observed  in  the  third  stage  of 
the  experiment.  Nonetheless,  the  uncertainty  in  the  recorded 
temperature  should  be  small,  compared  to  the  differential 
between  the  receiver's  external  case  and  the  oscillator. 

The  difference  in  patterns  shown  by  the  receivers, 
might  well  be  an  effect  of  an  uneven  heat  exchange.  Receiver 
2  had  to  be  placed  on  top  of  the  receiver  1,  because  not 
enough  room  was  available.  When  they  were  removed  from  the 
refrigerator,  receiver  1  was  warmer,  to  the  touch,  than 
receiver  2.  There  was  another  problem  with  this  experiment: 
two  hours  in  the  refrigerator  was  not  enough  for  the  receivers 
to  reach  a  steady  state. 

A  net  effect  of  the  temperature  on  the  frequency  is 
evident  now.  In  an  attempt  to  quantify  the  effect  of  the 
temperature,  a  new  experiment  was  conducted  on  August  11th. 
Both  receivers  and  temperature  sensors  were  left  a  whole  night 
inside  the  refrigerator.  This  time,  receiver  2  was  below  the 
receiver  1,  with  about  an  inch  spacing  between  them. 

After  one  and  a  half  hours  of  logging  data,  all  the 
hardware  was  removed  from  the  refrigerator.  They  were  left  at 
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Figure  8.  Thermal  shock  experiment,  08.08.91 
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ambient  temperature  for  two  and  a  half  more  hours.  Figure  9  is 
a  plot  of  the  Doppler  residuals,  from  a  linear  model,  and 
recorded  temperature.  In  this  case  the  steady  state  was 
reached,  both  inside  and  outside  the  refrigerator.  Receiver  2 
shows  a  slightly  smaller  variation  in  frequency.  This  can  be 
due  to  either  different  thermal  constant,  or  different  rate  of 
heat  exchange  inside  the  refrigerator.  The  crystal  oscillators 
are  very  likely  to  have  different  thermal  constants;  they  are 
not  thought  to  be  subject  to  very  stringent  specifications  due 
to  the  low  price  of  the  receiver. 

The  measured  variation  in  frequency  and  temperature 
for  receiver  2  leads  to  a  thermal  constant  of  16.3  m/s/degree 
in  the  second  test,  while  a  value  of  14.9  m/s/degree  was 
obtained  in  the  first  experiment,  for  the  same  receiver. 
Steady  state  was  not  reached  then.  These  results  are  in  good 
agreement  if  some  adjustment  is  allowed  accounting  for  this. 
Data  for  receiver  1  in  the  first  experiment  cannot  be 
considered  valid.  A  thermal  constant  of  18.9  m/s/degree  was 
obtained  in  the  second  experiment  for  this  receiver. 

As  a  conclusion,  the  thermal  dependence  of  the 
receiver  clock  can  be  estimated  to  be  somewhere  between  15  and 
20  meters  per  second  per  degree  centigrade.  This  is,  between 
0.5  and  0.7  hertz  per  degree  centigrade,  relative  to  the 
nominal  frequency,  10.23  Mhz. 
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VI.  CHARACTERI2AT10N  OF  THE  CLOCK  NOISE 


The  accuracy  of  a  clock  can  be  expressed  as  the 
combination  of  a  deterministic  and  a  random  component.  Stein 
[Ref.  12 :p.  207]  and  Allan  [Ref.  13 :p.  3]  state  that  precision 
oscillators  can  be  described  by 

x(t)  =Xo+yot+lDt^+6(t)  (63) 


where  x(t)  is  the  offset  of  the  clock,  relative  to  the 
absolute  time.  The  first  three  terms  in  eg. (63)  are  the 
deterministic  component  and  €(t)  is  the  random  component.  X(, 
represents  the  initial  clock  offset,  and  is  termed 
synchronization  error,  represents  the  initial  fractional 
frequency  offset  Af/fg  relative  to  the  nominal  frequency.  It 
is  known  as  syntonization  error.  D  is  the  variation  in  the 
fractional  frequency  error  and  is  known  as  frequency  drift  or 
aging. 

The  degree  to  which  this  second  order  polynomial 
successfully  models  the  deterministic  component  of  the 
receiver  oscillator  is  a  function  of  the  specific  oscillator 
type  and  quality.  For  quartz  crystals,  y^,  is  affected  by  the 
environmental  conditions  (see  V.C.2)  and  D  may  not  be  a 
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constant.  For  the  case  at  hand,  unless  the  environmental 
effects  are  completely  isolated,  it  is  not  possible  to  model 
the  deterministic  behavior  of  the  clock.  Thus,  these  effects 
will  mask  a  great  deal  of  the  random  noise. 

The  following  is  a  theoretical  discussion  of  the 
characterization  of  the  random  component  of  the  clock  offset. 
This  background  makes  it  possible  to  address  the  problem  of 
characterizing  the  noise  in  the  MX4200  receiver. 


A.  MEASURING  THE  FREQUENCY  STABILITY  IN  THE  TIME  DOMAIN 

Assume  an  oscillator  generates  a  voltage  signal: 


V(t)=Vo  sin*(t) 
=Vo  sin(2?rfot+</>(t) ) 


(64) 


Here  fg  is  the  nominal  frequency  and  <t>(t)  is  the  residual 
phase.  This  residual  phase  accumulates  the  effect  of  the 
instantaneous  deviation  of  the  actual  frequency  f (t)  from  fp. 
Because  measuring  equipment  of  infinite  bandwidth  is  not 
physically  feasible,  f (t)  is  not  a  measurable  quantity. 

While  the  instanatneous  frequency  is  not  measurable,  it  is 
the  starting  point  of  the  timing  analysis.  It  is  defined  as 


f,  1  d*{t) 
2n  dt 


(65) 
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Because  the  deviation  from  fg  can  be  very  small,  a 


fractional  frequency  offset  has  been  defined  as 


y(t)  = 


1  d0(t) 

2ir/o  dt 


(66) 


It  has  the  added  advantage  that  it  makes  the  analysis  of  the 
frequency  fluctuations  independent  of  any  scaling  factor 
between  the  oscillator  and  the  particular  circuit  where  the 
measurements  are  being  taken. 

Integrating  eg.  (66)  yields  the  instantaneous  clock  offset: 

I 

Jf(t)  =  fy(T)dT  =  .^iill^l£I  (67) 

i  27rfo 

When  the  signal  is  sampled  at  discrete  intervals,  an 
average  or  mean  fractional  frequency  offset  at  epoch  t*  can  be 
defined.  This  is, 

y*=Y  I 


where  r  is  the  sampling  interval. 
Then,  from  eg. (67) 


y*= 


x(t*+T)-x(tt) 

T 


(69) 
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1.  Th«  Allan  variance 

Certain  noise  processes  in  oscillators  have  a  standard 
deviation  that  does  not  converge  as  N-mo,  i.e.,  the  fractional 
frequency  offset  grows  unbounded  as  the  size  of  the  data  set 
increases  [Ref.  12:p.  201].  Then,  the  standard  deviation  is 
not  well-suited  for  long  term  analysis.  This  severely  limits 
its  applicability  for  characterizing  a  clock  behavior.  The 
long  term  stability  is  in  fact  what  determines  the  clock 
performance  as  a  time  keeper. 

To  deal  with  this  problem,  the  IEEE  recommended  the 
two-sample,  no  dead  time,  variance  of  the  mean  frequency 
offset,  as  a  better  defined  descriptor  [Ref.  13:p.  4].  It  is 
known  as  the  Allan  variance; 

Oy^iT)  =  <a/(2,T,T)>=  <70) 


It  has  been  shown  to  be  well  defined  for  all  the  major  noise 
processes  found  in  an  oscillator  [Ref.  12:p.  201]. 

Equations  (69)  and  (70)  allow  the  computation  of  the 
Allan  variance,  based  on  measured  time  offsets: 

aMT)=-L  <[x(t+2T)  -2x(t*T)  ^x{t)  ]2>  (71) 

2t^ 

Here  the  averaging  extends  to  an  infinite  time.  For  N  discrete 
time  readings,  a  good  estimate  of  the  Allan  variance  is: 
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1 


(72) 


N-2 

2(N-2)T^i 


53  [Jf(  V2T)  -2X(  Vt)  +X^(tJ  y 


Each  pair  xftj,  defines  a  single  average  fractional 
frequency  offset  in  the  interval  (t^,  t^+r)  . 

Because  no  dead  time  exists  between  measurements,  it 
is  possible  to  define  different  sampling  intervals  to  be 
integer  multiples  of  an  initial  r^.  Therefore 


dMmr)  = 


S-7m 


2  (N-2jn)m^T(,^ 


53  Cx(t*+2mTo) -2x(tt+JnTo) +x(tj  (73) 


gives  a  good  estimate  of  the  Allan  variance  for  a  sampling 
interval  mr^. 

2.  Confidence  intervals 

Consider  a  series  of  time  offset  measurements  x, 
spaced  t  in  time.  Two  consecutive  values  x*,  x*+;  yield  a  single 
mean  fractional  frequency  offset  y*.  x*+;  and  x^+j  define  y*+,. 
Thus,  those  three  consecutive  values  give  a  single  sample 
Allan  variance. 

For  the  common  oscillators,  the  first  difference  of 
the  frequency  is  a  normally  distributed  variable  with  zero 
mean  [Ref.  12;p.  210].  The  square  of  a  normally  distributed 
variable  follows  a  Chi-squared  distribution  with  one  degree  of 
freedom . 
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Pigur«  9.-  Different  single  sample  variances,  s,,  S2,_S3,  from 
the  same  data,  k  is  the  clock  reading  niimber,  and  y  is  the 
mean  fractional  frequency  error  corresponding  to  consecutive 
clock  readings 

Consider  now  the  five  time  offset  values  shown  in 
Figure  9.  From  readings  l  to  3  we  can  get  a  single  sample 
Allan  variance  from  readings  3  to  5  a  new  sample  Allan 
variance  s/(t)  is  obtained.  In  averaging  them,  a  better 
estimate  of  the  true  Allan  variance  is  obtained 

5/(T)=i(s,2+S22)  (74) 


which  has  a  tighter  confidence  interval.  Now,  9y{t)  ,  as  a  sum 
of  two  squared,  independent,  normally  distributed  variables, 
follows  a  Chi-squared  distribution  with  two  degrees  of 
freedom. 
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For  a  better  exploitation  of  the  data,  a  new  sample 
variance  s/ct)  can  be  obtained  as: 


SiHy)=^(y3~y2)^ 


(75) 


But  this  sample  variance  is  not  independent  of  the  previously 
computed  ones.  A  better  estimate  of  the  true  Allan  variance  is 
obtained  by  averaging  all  three.  But  this  estimate  does  not 
follow  a  Chi-squared  distribution  with  three  degrees  of 
freedom. 


The 

actual 

number 

of  degrees 

of 

freedom  is  a 

complicated 

matter. 

beyond 

the  scope 

of 

this  work.  It 

basically  depends  on  the  underlying  noise  type,  and  can  take 
fractional  values.  Barnes  [Ref.  14:p.  103  gives  expressions  to 
compute  the  degrees  of  freedom  as  a  function  of  the  number  of 
samples,  and  the  sampling  interval.  Those  expressions  are 
reasonably  well  suited  for  the  five  major  noise  types  found  in 
an  oscillator. 

From  the  point  estimate  computed  above,  an  interval 
estimate  of  the  Allan  variance  can  be  obtained  as  follows: 


OyUr) 


< 


».i- 


1 


at  a  confidence  level  a,  for  v  degrees  of  freedom. 


(76) 
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As  can  be  seen,  the  computation  of  this  interval 
estimate  is  somewhat  complicated.  First,  it  involves 
determining  the  number  of  degrees  of  freedom.  In  turn,  it 
implies  some  assumptions  about  the  kind  of  noise  present  in 
the  data,  for  the  given  sampling  interval.  In  addition,  the 
number  of  degrees  of  freedom  is  only  empirically  obtainable 
for  some  of  the  stochastic  processes  in  an  oscillator  [Ref. 
14 :p.  7].  Finally,  the  Chi-squared  distribution  for  a  non 
integer  number  of  degrees  of  freedom  has  to  be  used. 

A  simpler  approach  is  given  by  Stein  [Ref.  12:p.  213], 
citing  Lesage  and  Audoin  [Ref.  15]  and  Yoshimura  [Ref.  16].  It 
is  based  on  computing  the  variance  in  the  estimate  of  the 
Allan  variance. 

For  N  single  time  offset  values,  N-2  sample  Allan 
variances  are  obtained.  Averaging  them  all,  the  resulting 
estimate  dy('r)  will  deviate  AfW;  from  the  true  Allan  variance: 

a,^(T)=(l+A(N))a/  (77) 

Here  ^(N)  has  mean  zero.  For  N  greater  than  10,  the  variance 
of  A  is  approximately  given  by 

a^{A)«5:  (78) 
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where  C^,  again,  depends  on  the  kind  of  noise  in  the  process. 
Thus,  the  true  Allan  variance,  at  a  67%  confidence  level,  can 
be  approximated  by 


&y(r) (1 


■15 


)  <  0(T)  < 


SAT)  (1  + 


(79) 


Table  4  shows  the  value  of  for  the  five  major  noise 
processes  in  an  oscillator. 


TABLE  4.  VALUES  OF  FOR  THE  MAJOR  NOISE  TYPES 


Noise  Type 

a 

White  Phase 

2 

3.88 

Flicker  phase 

1 

3.88 

White  frequency 

0 

2.99 

Flicker  frequency 

-1 

2.31 

Random  walk  frequency 

-2 

2.25 

B.  STOCHASTIC  PROCESSES  IM  AN  OSCILLATOR 

The  stochastic  clock  noise  is  often  represented,  in  the 
frequency  domain,  as 


CO 


(80) 


Where  f  is  the  Fourier  frequency.  Sy(f)  is  the  one  sided  power 
spectral  density  of  the  phase  fluctuations,  and  h„  is  the 
amplitude  of  the  a-th  term. 
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Five  major  stochastic  processes  have  been  described  in 
precision  oscillators  [Ref.  13:p.  4] [Ref.  12;p.  204].  They 
correspond  to  the  terms  -2  to  2  and  named  as  follows: 

•  Random  walk  noise,  frequency  modulation  (FM) .  a=-2. 

•  Flicker  noise,  FM.  a=-l. 

•  White  noise,  FM.  a=0. 

•  Flicker  noise,  phase  modulation  (PM) .  a=l. 

•  White  noise,  PM.  a=2. 

Those  five  stochastic  processes  also  can  be  described  in 
the  time  domain.  A  plot  of  the  square  root  of  the  Allan 
variance,  in  a  log/ log  scale,  typically  looks  as  shown  in 
Figure  10.  The  sharp  corners  between  segments  is  a 
simplification  for  discussion  purposes. 

A  different  kind  of  noise  dominates  in  each  of  these 
regions.  The  slope  of  the  corresponding  segment,  n  is  given 
below: 

•  I.-  White  PM  and  Flicker  PM,  n»-l. 

•  II.-  White  PM,  n=-ll2. 

•  III.-  Flicker  FM,  n=0. 

•  IV.-  Random  walk  FM,  n=+l/2. 

This  pictorial  description  is  very  useful  in  characterizing 
the  noise  in  a  frequency  standard. 

Region  III  shows  no  dependence  on  the  sampling  interval. 
It  is  known  as  the  Flicker  floor.  Its  physical  cause  is  not 
fully  understood  yet.  It  represents  the  limit  in  the  frequency 
stability  of  a  clock,  i.e.,  its  accuracy.  For  a  time  standard. 
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If 


Figure  lO.-  A  typical  Allan  variance  plot 

it  is  desirable  to  have  the  flicker  floor  at  the  longest 
sampling  interval  possible.  The  best  measured  frequency 
stability  reported  by  Hellwig  [Ref.  17;p.  566]  for  passive 
hydrogen  maser  resonators  locates  the  Flicker  floor  at  a 
sampling  interval  longer  than  10*  seconds. 

For  crystal  oscillators,  that  Flicker  floor  can  be  located 
somewhere  between  10'^  and  10^  seconds.  Thus,  they  are  not  well 
suited  for  a  long  term  time  keeping.  For  short  term  intervals, 
they  can  show  a  better  stability  than  most  atomic  resonators, 
though.  That  is  why  both  kinds  of  oscillators  are  used 
together  to  form  a  time-frequency  standard. 

Region  IV  is  extremely  sensitive  to  the  frequency  drift  D 
in  eq. (63)  .  The  clock  offset  grows  at  a  rate  of  the  square  of 
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the  sampling  interval  if  D  is  not  zero.  A  linear  frequency 
drift  would  yield  a  slope  m=1  in  region  IV.  As  a  deterministic 
effect,  it  should  be  removed  from  the  data  beforehand.  This  is 
not  an  easy  task  because  the  drift  is  highly  dependent  on 
environmental  effects,  and  its  effect  is  integrated  over  a 
very  long  period. 

C.  CEARACTERI2ATZON  OF  THE  NOISE  IN  THE  MZ4200  RECEIVER  CLOCK 

The  preceding  discussion  provides  the  theoretical 
background  for  the  analysis  of  the  noise  in  the  receiver 
clock.  The  common  approach  for  this  analysis  is  to  compute  the 
clock  offset  at  many  different  epochs,  spaced  t.  Eq. (63) 
models  this  offset  as  a  function  of  the  interval  following 
synchronization.  The  instantaneous  offset  x(t)  has  to  be 
computed  as  the  difference  between  readings  of  the  clock 
(under  investigation)  and  a  frequency  standard.  It  is  assumed 
that  the  frequency  standard  has  such  a  high  degree  of 
frequency  stability  that  the  effect  of  its  fluctuations  on 
x(t)  is  negligible. 

The  GPS  SV's  provide  such  a  frequency  standard,  in 
comparison  to  a  crystal  oscillator.  The  readings  of  the 
receiving  clock  can  be  compared  with  this  standard  by 
subtracting  the  true  range  from  the  measured  range,  as  was 
explained  in  Chapter  V.  The  program  MXCKMO  already  corrects 
the  receiver  clock  for  the  major  deterministic  effects. 


90 


including  its  own  model.  Then,  the  output  residuals  are  in 
fact  €  (t)  in  eq.  (63)  . 

Here,  the  aforementioned  assumption  about  the  standard's 
frequency  stability  is  legitimate.  The  noise  level  in  the  SV 
clock  is  several  orders  of  magnitude  below  the  noise  in  the 
receiver  clock.  Several  systematic  factors  in  the  whole 
measurement  process  may  introduce  noise  in  the  residuals, 
though.  Those  in  the  higher  frequency  range  would  mask  the 
actual  noise  in  the  clock  at  short  sampling  intervals.  In 
particular: 

•  Pseudorange  noise 

•  Changes  in  the  selected  ephemeris 

•  SA  activation 

•  Loss  of  lock  on  the  SV 

need  to  be  carefully  watched.  Other  long  term  deterministic 
effects  were  discussed  in  Chapter  V.  The  effect  of  the 
temperature  on  the  long  term  noise  will  be  seen  here,  too. 

1.  Analysis  tools 

a.  Program  MX ALLY 

The  main  analysis  tool  here  is  the  program  MXALLV, 
listed  in  Appendix  I.  It  computes  the  Allan  variance  of  the 
range  residuals  following  eq. (73) .  Programs  MERGE  and  MXFXRS 
act  as  interfaces  between  the  MXCKMO  output  and  MXALLV. 
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b.  Program  MERGE 

This  program  generates  a  single  file  with  residuals 
from  different  satellites.  The  maximum  time  a  SV  can  be 
tracked  is  about  five  hours.  This  is  the  time  the  SV  is  above 
the  observer's  horizon.  With  this  program,  a  single  residuals 
file  can  be  obtained  for  a  longer  observation  session.  It 
accepts  residuals  if  they  are  not  flagged  as  anomalous.  Once 
it  gets  to  the  end  of  the  file,  or  when  a  flag  is  detected, 
the  program  prompts  a  new  input  file.  Then,  it  fetches  the  new 
file  until  the  last  valid  epoch  in  the  previous  file  is  found. 
The  difference  between  residuals  of  both  SV's  will  be  applied 
to  all  the  residuals  in  the  new  file.  This  allows  to  get  a 
continuous  sequence  of  residuals  from  the  beginning  to  the  end 
of  the  session.  Otherwise,  discontinuities  in  that  sequence 
would  degrade  the  Allan  variance  at  short  sampling  intervals. 

c.  Program  MXFXRS 

Gaps  may  exist  in  the  residuals  file  that  MXCKMO 
outputs.  For  the  Allan  variance  to  be  properly  computed,  a 
continuous  sequence  of  samples  has  to  be  taken.  If  the  number 
of  missing  epochs  is  small,  MXFXRS  synthesizes  the  missing 
residuals  by  interpolating  between  the  closest  records.  This 
approach  is  not  valid  if  more  than  1%  of  the  data  is  missing, 
because  it  artificially  decreases  the  noise  level. 
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d.  Program  CONFIN 

This  program  makes  it  possible  to  go  from  the  point 
estimate  computed  by  MXALLV  to  the  interval  estimate  given  in 
eq. (76) .  The  user  selects  a  confidence  level,  and  the  point 
estimate  Allan  variance  is  then  displayed.  The  decision  about 
the  boundaries  between  noise  regions,  as  depicted  in  Figure 
10,  is  left  to  the  user.  The  number  of  degrees  of  freedom  is 
computed  accordingly,  using  the  equations  in  Barnes  [Ref. 
14;p.  10].  The  parameter  a  in  eq. (76)  is  being  shifted  until 
the  probabilities 


and 


Then,  the  interval  extremis  are  found  as 


Oy(T)=dAr) 


a,<x/2 


(82) 


(83) 


^^nl-a/2 


(84) 
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2.  Short  torn  analysis 

The  effect  of  the  temperature  on  the  clock  noise 
should  not  be  a  major  problem  for  short  periods.  An  example  is 
examined  here  to  confirm  this.  An  observation  session  was  held 
on  June  18th,  1991  (day  169),  under  the  following 

circumstances : 

•  Receiver  1. 

•  Time  recovery  enabled. 

•  Data  logged  at  1  s.  rate,  during  3.3  hours. 

•  The  receiver  was  inside  a  faixly  well  insulated  room. 

•  SA  was  not  activated. 

A  quadratic  clock  offset  model  was  fitted  to  the  data. 
That  meant  that  a  constant  drift  in  the  clock  frequency  was 
assumed.  If  the  clock  frequency  did  in  fact  have  a  constant 
drift,  the  resulting  Allan  variance  accurately  describes  the 
stochastic  noise.  No  temperature  record  was  kept  in  that 
session.  A  constant  drift  in  the  frequency  might  have  been  due 
to  either  a  constant  rate  of  change  in  the  temperature,  or 
aging.  It  does  not  seem  reasonable  to  expect  aging  effects  for 
such  a  short  period.  If  there  were  any,  they  would  be  masked 
by  temperature  effects.  To  assume  a  linear  change  in 
temperature  over  this  period  seems  to  be  an  over 
simplification. 

The  resulting  Allan  variance  is  shown  in  Figure  11. 
The  general  shape  of  the  curve  matches  fairly  well  with  what 
should  be  expected.  The  Flicker  floor  is  located  at  about  one 
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Rec  nl  day  l<-9,  residuals  off  a  quadratic  model 
<90  V.  confidence  interval) 


Figure  11.-  The  short  term  Allan  variance  in  receiver  1 

second.  In  the  range  1-20  seconds,  both  Flicker  FM  and  Random 
Walk  FM  noises  are  present.  The  last  dominates  at  sampling 
intervals  longer  than  20  seconds.  There  are  many  crystal 
oscillators  that  show  a  Flicker  floor  all  the  way  down  to  the 
one  second  sampling  interval.  Shorter  samples  could  not  be 
taken  because  one  second  is  the  interval  between  measurement 
epochs . 

Another  observation  session  was  held  on  May  30th,  1991 
(day  150) ,  as  follows: 

•  Receiver  2 

•  Time  recovery  enabled 

•  Data  logged  at  1  second  rate,  during  40  minutes 
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•  SA  was  not  activated 


•  Environmental  conditions  expected  to  be  as  stated  before 
for  receiver  1. 

A  quadratic  model  was  fitted  to  the  residuals  from  SV 
6.  In  this  case,  the  observation  session  was  much  shorter, 
although  Allan  variance  up  to  one  thousand  seconds  was 
obtainable.  The  confidence  interval  is  wider,  too.  They  are 
plotted  in  the  upper  part  of  Figure  12. 

In  comparing  the  results  for  both  receivers,  a  smaller 
level  of  noise  was  noticed  for  the  number  2.  There  is  also  a 
qualitative  difference.  The  Flicker  floor  for  receiver  2  is 
somewhere  between  l  and  4  seconds  sampling  interval  and  white 
FM  noise  seems  to  dominate  at  the  measurement  interval,  one 
second.  The  noise  level  at  this  sampling  interval,  for 
receiver  2,  is  about  2/3  the  one  in  receiver  1.  It  is  less 
than  1/2  for  sampling  intervals  longer  than  4  seconds. 

To  confirm  these  results,  data  from  a  previous 
observation  session  were  analyzed.  They  were  gathered  the  day 
before,  under  similar  circumstances,  with  receiver  2.  The 
Allan  variance  of  the  phase  range  residuals  to  SV  6  is  plotted 
at  the  bottom  of  Figure  12.  A  quadratic  clock  model  was 
removed ,  too . 

In  this  case,  a  non-constant  linear  drift  seems  to  be 
present  in  the  clock.  It  could  be  due  to  temperature  changes. 
The  slope  in  the  tail  of  the  curve  is  greater  than  1/2.  That 
means  that  deterministic  effects  are  masking  the  long  term 
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random  walk  FM  noise.  Nevertheless,  those  effects  should  be 
negligible  at  sampling  intervals  of  a  hundred  seconds  or  less. 
In  fact,  the  upper  and  lower  parts  of  Figure  12  show  a  very 
good  agreement  in  the  short  term  noise  levels. 

In  contrast,  noise  for  session  150,  receiver  1, 
(Figure  13)  shows  the  same  characteristics  as  we  found  before. 
The  set  up  was  the  same  as  for  receiver  2.  The  results  are 
very  consistent,  no  matter  what  the  length  of  the  session  was 
(40  minutes  versus  3  hours  and  20  minutes) .  In  conclusion,  a 
higher  noise  level  in  receiver  1,  on  a  short  term  basis,  has 
been  measured. 


R«c  til  1S0.  r«sidu*Is  off  •  qu«dr*iie  model 

<  y.  eorifidence  interw«il> 


Figure  13.-  Allan  variance  for  receiver  1,  day  150. 
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3.  High  frequency  systematic  effects 

Four  major  systematic  effects  were  mentioned  above 
that  should  be  carefully  watched  in  getting  the  Allan 
variance. 

a.  Pseudorange  noise 

The  residuals  that  the  clock  model  program  outputs 
are  computed  using  phase  derived  ranges.  At  short  intervals, 
higher  lev->ls  of  noise  should  be  expected  using  pseudoranges, 
because  of  the  nature  of  the  measurement  process. 

A  slight  modification  in  the  MXCKMO  program  was 
done  to  get  the  residuals  from  the  pseudoranges,  instead  of 
phase  ranges.  The  resulting  Allan  variance  is  plotted  in 
Figure  14  as  (x)  .  The  corresponding  Allan  variance  of  the 
phase  range  residuals  is  plotted  as  (+) .  The  noise  here  is 
dominated  by  processing  noise  in  the  correlator.  Its  level  is 
more  than  one  order  of  magnitude  at  the  shorter  sampling 
times.  Both  of  them  converge  for  sampling  intervals  larger 
than  100  s.  This  would  be  a  good  time  constant  for  the  CODE 
variable  in  eg. (42)  to  be  filtered,  before  being  input  into 
the  navigator. 

b.  Ephemeris  changes 

This  is  not  a  trivial  question,  sometimes  not 
properly  accounted  for.  The  parameters  needed  for  the  true 
range  computation  are  contained  in  the  broadcast  ephemeris. 
Different  ephemeris  sets  are  chosen  as  a  function  of  the  time 
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Ree  t!  1  4ai<i  )<.9,  noise  In 
Pseudorenge  <x>  end  pitesc  renge  <  +  > 


Figure  14.-  Processing  noise  in  the  pseudorange,  receiver  1 

difference  between  the  ephemeris  epoch  and  the  measurement 
time.  The  computed  range  to  the  SV,  then,  shows 
discontinuities,  whenever  the  program  switches  from  one 
ephemeris  set  to  another. 

This  has  a  noticeable  effect  on  the  high  frequency 
noise,  if  not  properly  filtered.  The  program  MXCKMO  biases  the 
OMC's  by  those  discontinuities,  so  the  resulting  true  range  is 
smoothed.  An  estimate  of  the  effect  of  leaving  these  jumps  in 
follows. 

Assume  a  discontinuity  S  in  the  OMC  exists  between 
epochs  j+1  and  j+2,  as  shown  in  Figure  15.  Two  sample 
variances  will  be  affected:  (j,  j+2)  and  {j+1,  j+3) .  The  OMC 
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Figur*  IS.-  Effect  of  an  ephemeris  change  on  the  'observed 
minus  computed'  (OMC)  range. 

is  represented  by  x  before  the  discontinuity,  x'  after  it. 
Following  eg.  (72)  ,  let  a^,'^  be  the  Allan  variance  corresponding 
to  these  OMC's: 


j-i 


2To  (W-2)  t.| 

+  (x^.^j+5-2  (x^.,2+5)  +x^„)2 

^-2 

+  E  (^f2-"«-2(x*.,  +  5) +x*+5)2] 

*•>♦2 


(85) 


All  the  5's  in  the  last  term  cancel.  Gathering  terms  in  x, 
eg. (85)  becomes 
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(86) 


+  2  «2+2  5  (X,>j-X>*2 -■*}*!  *Xj)  ] 

Let  Oy  be  the  Allan  variance  corresponding  to  the  unbiased 
ONC's,  i.e.,  to  the  x's  before  and  after  the  discontinuity. 
The  first  term  in  eg.  (86)  is  in  fact  this  value.  The  last 
parenthesis  in  eg. (86)  must  be  small,  compared  to  S.  Values  of 
S  larger  than  ten  meters  have  been  observed.  The  contribution 
of  the  last  term,  thus,  is  negligible. 

Eg. (86) ,  then,  becomes 


I  2  i 

y 


(W-2)ro^ 


(87) 


for  a  single  discontinuity.  If  M  discontinuities  exist  in  the 
data,  each  one  contributes  a  new  term.  It  should  be  pointed 
out  that  changes  in  ephemeris  can  happen  once  in  an  hour  for 
a  given  SV.  Thus, 


a'  ^»o  ^ 
y  y 


(W-2)To 


(88) 


Their  effect,  at  a  sampling  interval  rnTj,  will  be 
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(89) 


M 


<7'y  (fflTo)«a^(jnTo)+. 


although  as  m  grows,  the  last  term  in  eg.  (86)  may  become 
comparable  to  6^, 

Figure  16  shows  the  effect  of  changing  ephemeris  on 
the  Allan  variance,  at  short  time  intervals.  The  points  marked 
(x)  correspond  to  the  Allan  variance  of  the  data  if  those 
jumps  are  not  removed.  This  is,  in  eg. (89) .  Points  marked 
as  (+)  correspond  to  a^,  i.e.,  the  Allan  variance  of  the 
residuals,  once  filtered  by  discontinuities  in  OMC.  It  is,  in 
fact,  the  same  Allan  variance  shown  in  Figure  11.  Residuals 
from  SV  20  were  used  there.  A  change  in  ephemeris  occurred  at 
epoch  262814.  Its  effect  on  the  computed  range  had  a  magnitude 
of  >7.14  meters.  Let  us  check  eg. (88)  against  the  values  shown 
in  Figure  16.  In  this  case,  N  was  12053,  and  M  was  1.  For  a 
one  second  interval 


6^ 

T^HN-2m) 


* - _ =4.7-  10-20 

(3-  10*)2-  12051 


(90) 


where  the  increment  in  computed  range  has  been  scaled  to  time. 
From  the  plots  in  Figure  16  we  get: 

<^'y(ro)  -ff,^(fo)  =(2.2-  10-‘®)2-(0.49'  10-‘®)2=4.59-  lO’^  (91) 
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Rec  til  day  li!.9 
Effect  of  erl->eineris  elianges 


Figure  16.-  The  effect  of  ephemeris  changes  in  the  short  term 
Allan  variance 

For  a  different  sampling  interval,  the  agreement  is  not  as 
good.  If  m  =  5 


_ (-7-14)^ _ 

(12053-2*  5)  *  5^*  (3*  10*)^ 


=0.19* 


(92) 


whereas  from  the  plot 

a'y^(5) -Oy^(5)  =0.26*  10-^  (93) 
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As  the  sampling  interval  increases,  the  last  term  in  eg. (86) 
may  effectively  increase,  compared  to  S^.  It  produces  slightly 
different  results. 

c.  Temporary  loss  of  lock  and  cycle  slips 

A  tracked  SV  can  be  temporarily  lost.  At 
reacquisition,  a  different  bias  is  introduced  in  the  phase- 
derived  range.  This  bias  is  on  the  order  of  the  noise  in  the 
pseudorange,  i.e.,  several  meters.  Its  effect  on  the  Allan 
variance  is  expected  to  be  of  the  same  order  as  the  changes  in 
ephemeris.  Nothing  can  be  done  to  correct  for  the  new  bias. 
But  the  data  should  be  carefully  watched  to  avoid  computations 
on  residuals  from  a  SV  that  was  temporarily  lost.  The  program 
MXCKMO  can  detect  spikes  in  the  OMC  larger  than  about  a  meter. 
Thus,  a  reasonably  quantity  of  them  are  flagged  as  bad 
records.  The  Allan  variance  routine  will  not  accept  a 
residuals  file  with  those  kinds  of  records  in  it. 

This  problem  is  especially  harmful  when  the  logging 
interval  is  longer  than  one  second.  Otherwise,  the  lock  loss 
is  easily  detected  if  measurements  to  that  SV  are  missing  for 
a  few  seconds. 

Cycle  slips  have  not  been  detected,  and  it  appears 
to  be  very  unlikely.  Their  size  is  in  the  order  of  decimeters, 
so  it  would  be  hard  to  detect  them.  Their  effect  on  the  Allan 
variance  would  not  be  much  of  a  problem,  because  they  are  so 
infrequent. 
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d.  Selective  Availability  activation 


Not  too  much  is  knovm  about  the  way  SA  is 
implemented.  Higher  level  of  noise,  both  in  CA  code  and  in  LI 
frequency  should  be  expected.  Whether  that  level  is  above  or 
below  the  receiver  noise  could  not  be  checked. 

SA  was  activated  on  July  2nd,  1991  [Ref.  18]. 
Several  pieces  of  data  were  logged  after  that.  The  observed 
noise  in  the  phase  range  residuals,  as  depicted  by  the  Allan 
variance,  stayed  at  the  same  level  as  before. 

4.  Long  term  analysis 

Long  term  here  refers  to  observation  sessions  lasting 
more  than  a  complete  pass  of  a  satellite.  The  analysis  has  to 
be  performed  on  a  sequence  of  residuals  of  more  than  one  SV. 
Residuals  from  different  SV's  are  biased  to  those  of  the 
earliest  one,  as  explained  in  Vl.C.l.b.  So,  the  entire 
sequence  is  internally  consistent. 

Plenty  of  data  were  gathered  in  long  term  sessions. 
Unfortunately,  the  effect  of  the  temperature  changes 
contaminates  most  of  those  data  sets.  A  few  more  tests  were 
done  attempting  to  thermally  isolate  the  receivers.  Sometimes 
the  receivers  were  taken  inside  a  refrigerator,  in  an  attempt 
to  keep  the  temperature  unchanged.  Comparisons  between  inside 
and  ambient  temperatures  were  done.  The  temperature  inside  the 
refrigerator  followed  the  main  ambient  trend,  by  a  factor  of 
one  half.  In  addition,  acceleration  effects,  due  to  the 
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compressor  cycle,  added  an  extra  complication.  No  significant 
result  was  obtained  by  this  procedure. 

A  second  approach  was  followed.  Receivers  were  left  at 
ambient  temperature,  in  a  closed  room,  with  an  external  heat 
source.  A  temperature  record  was  used  to  estimate  a  period  of 
constant  rate  of  change  on  temperature.  A  linear  drift  in  the 
frequency  should  be  expected  then.  This  is,  D  in  eg. (63) 
should  be 

D=D(t)=i?o+^t  (94) 


so  eg. (63)  becomes 

x(t)  =Xo+yot  +  .iDot2+lD,t^+e(t)  (95) 


Thus,  a  cubic  model,  rather  than  quadratic,  would  be 
needed  to  remove  this  deterministic  effect.  This  was  done  on 
data  from  day  221  (August  10th,  1991) .  Receiver  2  was  logging 
data  for  24  hours,  at  a  15  second  rate.  A  linear  model  was 
first  computed  for  the  whole  session.  A  plot  of  the  Doppler 
residuals  from  that  model  is  shown  in  Figure  17.  From  that 
plot,  a  linear  frequency  drift  can  be  estimated  between 
between  00. oo  PST  and  09.00  PST,  August  10th.  It  matched  the 
temperature  records. 

Data  were  then  processed  between  those  epochs.  This 
time,  a  cubic  model  was  generated.  The  resulting  Allan 


107 


Receiver  a2,  Doppler  resitfuels  off  e  linear  model 


Figure  17.-  Doppler  residuals  off  a  linear  model.  Receiver  2, 
day  221 

variance  is  shown  in  Figure  18.  The  slope  in  the  curve  (1/2) 
corresponds  to  the  expected  random  walk  FM  noise.  Systematic 
effects  seem  to  have  been  removed.  The  agreement  with  plots  in 
Figure  12  in  the  overlap  region  is  good. 

A  comparable  test  for  receiver  1  was  not  possible. 
Earlier  data  were  reviewed.  A  linear  fit  to  data  from  day  117 
showed  an  almost  constant  frequency  drift  between  20.00  PST, 
April  27th,  and  07.00  PST  April  28th  (Figure  19).  No 
temperature  data  were  available  then.  Data  were  reprocessed 
for  that  interval,  and  a  quadratic  clock  offset  model 
generated.  The  Allan  variance  of  the  residuals,  from  that 
model,  is  shown  in  Figure  20.  Again,  a  1/2  slope  means  that 
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Rac  a2  dav  221,  rasitfuals  off  a  euwie  model 
<90  X  confidence  interval  > 


Figure  18.-  Long  term  Allan  variance.  Receiver  2,  day  2il 

the  bulk  of  systematic  effects  has  been  removed.  It  shows  good 
agreement  with  those  in  Figures  11  and  13. 

The  reader  is  warned  about  the  validity  of  these 
results.  A  great  deal  of  empiricism  is  underlying  this  last 
procedure.  A  more  rigorous  treatment  of  the  temperature 
effects  is  needed.  In  particular,  an  environmental  chamber  is 
required  to  completely  isolate  the  noise  in  the  oscillator, 
from  those  induced  by  the  temperature.  In  addition,  analysis 
of  longer  data  sets  is  desirable.  Surveying  observation 
sessions  usually  last  for  some  three  hours.  The  noise  level  at 
that  sampling  interval  would  be  significant. 
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Receiver  til,  Doppler  residuals  off  e  linear  model 


Figure  19.-  Doppler  residuals  off  a  linear  model.  Receiver  1, 
day  117 
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Rec  al  dav  117.  residuals  off  a  quadratic  model 


<90  a  confidence  interval  > 


Figure  20.-  Long  term  Allan  variance.  Receiver  1,  day  117 
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VII.  CONCLUSIONS 


A  new  generation  of  inexpensive  GPS  receivers  is  emerging, 
which  utilize  temperature-compensated  crystal  oscillators.  The 
use  of  an  inexpensive  clock  implies  higher  rates  of  drift  and 
sensitivity  to  environmental  effects  than  customary  in  GPS 
applications.  The  guestion  of  the  clock  guality  in  one  of 
these  new  receivers,  the  Magnavox  4200,  has  been  addressed. 
The  characteristics  of  the  random  noise  in  the  measurements 
have  been  established  for  two  particular  receivers  by 
comparison  with  the  satellite  clocks.  The  effect  of  the 
temperature  on  the  oscillator  frequency  has  been  measured  as 
well . 

The  receiver  architecture  is  designed  to  minimize  the 
noise  in  the  pseudorange  by  means  of  a  Doppler  aided  code 
tracking  loop.  Doppler  frequency  wipe  off  occurs  before 
correlation.  By  integrating  this  Doppler  frequency,  a  phase 
derived  range  is  obtained.  The  time  delay  corresponding  to 
this  phase  derived  range  is  applied  to  the  code  replica 
generator.  Thus,  only  the  difference  between  pseudorange  and 
phase  range  has  to  be  tracked.  This  difference  has  a  small 
variation,  that  can  be  tracked  in  a  narrow  bandwidth. 

Measurements  are  simultaneously  taken  on  six  channels,  at 
receiver  time  epochs.  Clock  and  frequency  offset  are  easily 
obtained  as  a  part  of  the  navigation  solution,  using  an  eight- 
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state  Kalman  filter.  For  a  clock  with  such  a  high  variability, 
those  epochs  can  be  offset  several  milliseconds  from  the  GPS 
time  after  a  few  hours.  In  Time  Recovery  mode,  the  clock 
solution  is  used  to  dynamically  correct  the  internal  clock  on 
a  measurement  by  measurement  basis. 

Even  temperature  compensated  crystal  oscillators  are  very 
sensitive  to  temperature  changes,  at  the  levels  of  accuracy 
used  in  GPS.  The  effect  of  the  temperature  on  the  frequency 
offset  was  estimated.  Experiments  done  on  two  receivers, 
showed  that  the  transient  regime  lasted  for  two  hours.  Other 
tests  showed  that  the  thermal  constant  was  about  20  m/s/^C. 
Even  in  a  mild  environment,  this  may  easily  produce  clock 
offsets  of  several  miliseconds  in  just  a  few  hours. 

Temperature  effects  greatly  complicated  the  analysis  of 
the  stochastic  noise.  A  careful  process  was  followed  to 
isolate  deterministic  effects  from  the  measurements.  The 
stochastic  noise  in  two  receivers  was  obtained  for  sampling 
intervals  up  to  4000  sec.  They  show  very  similar 
characteristics.  The  Flicker  floor  is  at  about  the  measurement 
interval,  one  second.  The  Allan  variance  at  one  second 
sampling  interval  is  on  the  order  of  5x10'".  This  implies 
inter-measurement  random  errors  on  the  order  of  1.5  cm  due  to 
clock  offsets.  The  random  walk  frequency  noise  dominates  at 
intervals  of  4  seconds  or  more.  At  a  30  second  sampling  rate, 
the  Allan  variance  is  about  10'‘°,  implying  about  one  meter  of 
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clock  induced  random  error  in  range.  The  trends  imply  that 
over  a  three  hour  observation  session,  the  random  component  of 
the  range,  due  to  the  clock  noise,  would  grow  to  a  few 
kilometers.  The  time  tags  would  have  random  errors  of  tens  of 
microseconds . 

For  high  quality  geodetic  surveys,  this  has  an  important 
implication.  The  random  noise  can  accumulate  half  a  wavelength 
in  just  a  few  seconds,  making  integer  cycles  ambiguity 
resolution  by  clock  modelling  impossible.  For  all  intended 
purposes  the  clock  error  must  be  solved  for  on  a  measurement 
by  measurement  basis.  An  observation  session  would  require  a 
continuous  tracking  of  at  least  four  satellites,  with 
measurements  every  second. 
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APPENDIX  A.  AUTEOR'8  SOFTWARE 


1. -  MXCKMO  Computation  of  the  coefficients  of  a  polynomial, 

up  to  4th  degree,  that  models  the  receiver  clock 
offset.  Generates  phase  range  residuals  off  that 
model.  Linked  to: 

a. -  Subroutines  in  GAPP5  library:  ASKSTA,  ASKWTR,  GETEPH, 

LDBTAB,  RETRHO,  SETGPS,  STATL,  STATR,  STATXY,  WGS84. 

b. -  Subroutines  in  GENS  library:  WAIT. 

c. -  Functions  in  GAPP5  library:  ANLTRP,  ELVANG. 

d. -  Own  subroutine:  MXRDRl:  Data  interface  to  retrieve 

measurements  and  clock  shift.  Linked  to: 

1. -  Dr.  Clynch's  subroutines  in  MX42RTN  library: 

IM4SET,  M42SET,  NXMX42. 

2. -  Author's  generated,  following  [Ref.  14 :p.  31- 

39]:  LUBKSD,  LUDCMD,  MPROVD. 

2. -  MXFXRS  To  fill  in  missing  records  in  a  residuals  file; 

Halts  process  when  a  flagged  residual  is  found. 

3. -  MERGE  To  add  residuals  files  from  different  SV's  to  get 

a  continuous  sequence  of  residuals. 
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4. -  MXALLV  To  compute  point  estimate  the  Allan  variance  of 

the  residuals  from  a  model. 

5. -  MXTRND  To  remove  a  clock  offset  model  from  a  set  of  range 

measurements.  Linked  to: 

a. -  Subroutines  in  GAPP5  library:  ASKSTA,  ASKWTR,  GETEPH, 

LDBTAB,  RETRHO,  SETGPS. 

b. -  Functions  in  GAPP5  library:  ANLTRP,  ELVANG. 

c. -  Own  subroutine:  MXRDRl. 

6. -  CONFIN  To  compute  interval  estimate  Allan  variance,  from 

a  point  estimate  file.  Linked  to  the  following 
subroutines,  adapted  from  [Ref.  ll;p.  160-165]: 
GSER,  GAMMLN,  GCF  and  functions  GANMP,  GAMMQ. 
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APPENDIX  B.  ARL  SOFTWARE 


1. -  Executable  programs: 

a. -  DIFDAT2  To  get  differences  across  residuals. 

b. -  EPHSXJM  To  get  summary  of  an  ephemeris  file. 

c. -  GS2FIC  To  convert  an  ephemeris  file  in  FIC  format  to 

ASCII  format. 

2. -  Called  subroutines. 

a.-  GAPP5  library: 

1. -  ASKSTA  prompts  for  ground  station  input. 

2. -  ASKWTR  prompts  for  metorological  data. 

3. -  GETEPH  fetches  ephemeris  table  to  retrieve  the 

closest  set. 

4. ~  LDBTAB  Loads  array  of  ephemeris  in  a  file. 

5. -  RETRHO  computes  slant  range  to  the  SV. 

6. -  SETGPS  loads  GPS  parameters. 

7. -  STATL  Initializes  statistics  array. 

8. -  STATR  Inputs  new  value  to  statistics  array. 

9. *>  STATXY  Get  statistics  from  array  of  values. 

10. -  WGS84  Loads  WGS  84  reference  system  parameters. 

11. -  Function  ELVANG  computes  elevation  angle  of  a  SV. 

12. -  Function  ANLTRP  computes  Hopfield  tropospheric 

delay  model. 


119 


APPENDIX  C.  DR.  CLYNCH'S  SOFTWARE 


Executable  programs: 

a. -  GP:  General  Plot  package,  microcomputers  version. 

b. -  MX2EPH:  Extracts  ephemeris  from  a  data  file  and 

convert  them  into  a  FIC  binary  format. 

c. -  MXR  package:  data  logger  for  the  MX4200  receiver. 

d. -  MX42RTN  library;  contains: 

1. -  NXMX42  Reads  raw  data  files,  one  record  at  a  time, 

and  loads  values  in  arrays. 

2. '-  MXCRK2  Decompresses  block  2's  and  load  values  in 

arrays,  as  if  done  for  block  I's. 

3. -  IM4SET  Initialize  data-solution-time  block  arrays. 

4. -  M42SET  Gets  in  a  single  block  all  data,  solution 

and  time  records  belonging  to  the  same 
epoch . 
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APPEMDIX  D.  0ML8  ADJUSTMENT 


At  a  given  epoch  i  the  physical  clock  offset  is  first 
estimated  by  comparing  measured  pseudoranges  with  computed 
ranges.  These  observed-minus-computed  ranges,  scaled  to  time 
(e),  enter  in  an  Ordinary  Mean  Least  Squares  Adjustment 
(OMLS) .  The  mathematical  model  is: 


T  (t)  =a+b(t-to)  *£  (t-to)2+|  (t-to)’ 

=a+bAt+£At^+— 

2  6 


(96) 


X  •  G  •  / 

“(t)  =  H(At)  X  (97) 

Here  x  is  the  column  vector  of  the  polynomial  coefficients 
a,  b,  cf2,  d/6,  and 

H(At)  =  [l  At  At^  At’]  (98) 


is  the  design  matrix.  The  elements  At,  are  the  intervals 
between  measurement  and  reference  epoch,  at  epoch  i,  in  GPS 
time. 

The  parameters  are  first  set  to  0,  so  at  epoch  i 
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(99) 


Subsequent  iterations  lead  to  better  estimates  of  the 
parameters.  Once  the  convergence  criteria  is  met,  after  j 
iterations,  the  adjusted  coefficients  of  the  polynomial  will 
be: 
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(104) 


xi  =  jp/-‘  + 

The  whole  process  usually  converges  after  four  iterations. 
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